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Abstract 


This  study  compared  the  results  of  two  computer 
programs,  a  flux-difference-splitting  (FDS)  Godunov-based 
scheme  and  the  SCRAM jet  Hypersonic  Nozzle  (SCHNOZ) 
parabolized  Navier-Stokes  code  using  MacCormack's  method, 
applied  to  a  hypersonic  nozzle  flowfield.  Two  different 
nozzle  geometries  were  investigated  for  three  different  Mach 
numbers  along  a  typical  hypersonic  flight  trajectory.  A 
direct  comparison  between  the  SCHNOZ  and  FDS  programs  was 
made  by  numerically  solving  the  steady  Euler  equations  using 
a  frozen  flow  assumption  in  the  nozzle.  Significant  areas 
of  interest  in  comparison  of  the  code  results  were  the 
accuracy  in  capturing  the  flow  physics  and  the  required 
computational  time.  The  frozen  flow  SCHNOZ  program  is 
currently  6  to  10  times  more  efficient  in  terms  of 
computational  time  than  the  FDS  frozen  flow  program.  The 
SCHNOZ  and  FDS  codes  demonstrated  comparable  accuracy  in 
capturing  the  flow  physics  of  the  nozzle  flowfields 
considered.  The  implementation  of  the  viscous  terms  in  the 
SCHNOZ  code  proved  to  be  ineffectual  in  modeling  the  viscous 
effects  in  the  flowfield.  The  finite-rate  chemistry  effects 
were  important  for  the  nozzle  inlet  conditions  considered, 
as  the  SCHNOZ  finite-rate  chemistry  model  calculated  nozzle 
wall  thrusts  up  to  4%  greater  than  the  frozen  flow  model. 
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A  Comparison  of  Computational  Fluid  Dynamics  Computer 
Programs  for  Hypersonic  Propulsive  Nozzle  Flowfields 


I .  Introduction 


1 . 1  Background 

The  National  Aerospace  Plane  (NASP)  and  other  trans- 
atmospheric  vehicles  have  rejuvenated  interest  in  the 
hypersonic  flight  regime.  The  currently  proposed  propulsive 
system  in  the  hypersonic  regime  for  NASP-type  vehicles  is 
the  airframe-integrated  Supersonic  Combustion  RAM jet 
(SCRAM jet)  cycle.  Figure  1.1  shows  a  typical  hypersonic 
vehicle  with  an  airframe-integrated  SCRT^jet  nozzle.  The 
SCRAM jet  nozzle  contours  are  designed  to  generate  sufficient 
thrust  at  different  flight  conditions  given  a  certain 
vehicle  size,  weight,  and  fuel  onboard.  With  the  airframe- 
integrated  SCRAMjet,  this  in  turn  influences  the  engine 
size,  vehicle  size,  and  fuel  requirements  (13:50).  Numerous 
iterations  are  thus  required  in  the  optimization  of  the 
vehicle  and  engine  design. 

The  solution  to  the  hypersonic  propulsion  design 
problem  will  rely  heavily  upon  computational  fluid  dynamics 
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(CFD) .  This  reliance  on  CFD  has  been  driven  by  the 
decreasing  cost  of  computers  and  increasing  cost  of  wind 
tunnel  tests  coupled  with  a  limited  hypersonic  test  facility 
capability  (9:657).  The  hypersonic  propulsive  environment, 
with  associated  high  Mach  numbers  and  high  temperatures, 
reguires  very  accurate  CFD  simulations  (7:99). 

Additionally,  the  iteration  and  optimization  required  in  the 
design  process  of  a  hypersonic  propulsion  system  demands  an 
efficient  CFD  algorithm. 

1 . 2  Purpose 

Propulsive  nozzles  for  hypersonic  vehicles  represent  an 
extremely  demanding  test  for  CFD  codes.  The  flowfield  of  a 
propulsive  nozzle,  as  shown  in  Figure  1.2,  is  a  complicated 
structure  with  expansion  waves,  shock  waves,  contact 
surfaces,  and  the  interaction  among  all  three  and  surface 
boundaries  such  as  the  nozzle  wall  or  cowl  (6:1).  Several 
CFD  codes  have  been  developed  to  resolve  the  propulsive 
nozzle  flowfield  of  a  hypersonic  vehicle.  As  such,  the 
results  of  any  one  CFD  code  need  to  be  compared  with  those 
from  other  computer  codes. 

Doty  (6)  developed  a  flux-difference  splitting  (FDS) 
code  using  the  steady  Euler  equations  and  assuming  perfect 
gas  to  analyze  propulsive  nozzle  flowfields  in  an  effort  to 
optimize  nozzle  contours.  Doty's  FDS  code  was  later 
modified  by  Schieve  (13)  to  incorporate  a  calorically 
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imperfect  gas  thermodynamic  model.  A  SCRAMjet  Hypersonic 
Nozzle  (SCHNOZ)  parabolized  Navier-Stokes  (PNS)  computer 
code  was  developed  by  Science  Applications  International 
Corporation  (SAIC)  for  the  NASP  Program  (20) .  The  SCHNOZ 
code  includes  algorithms  for  both  viscous  flow  effects  and 
finite-rate  chemistry  reactions  in  numerically  solving  the 
PNS  eguations  for  the  flowfield.  The  purpose  of  this  study 
was  to  compare  results  from  the  modified  FDS  code  to  SCHNOZ 
code  results  to  determine  the  most  accurate  and  efficient 
code  in  the  analysis  of  a  hypersonic  propulsive  nozzle 
flowfield. 

1 . 3  Scope 

This  research  effort  sequentially  compared  the  accuracy 
and  efficiency  of  Doty's  FDS  code  with  the  thermodynamic 
model  improvement  of  Schieve  (13)  to  the  SCHNOZ  PNS  code 
(20)  applied  to  a  hypersonic  propulsive  nozzle  flowfield. 

The  nozzle  contour  in  the  study  was  a  two-dimensional, 
maximum  thrust,  planar  nozzle.  Included  in  the  study  were 
nozzle  contours  determined  from  Doty's  FDS  code  by  Herring 
to  be  the  optimum  nozzle  configuration  over  a  typical 
hypersonic  flight  trajectory  (8)  and  a  more  severe  expansion 
nozzle. 

The  flow  in  the  two-dimensional  hypersonic  nozzle  was 
assumed  to  be  steady,  compressible,  and  rotational.  The 
fluid  in  the  nozzle  was  treated  as  a  calorically  imperfect. 
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but  thermally  perfect,  gas  where  computations  of  the 
properties  downstream  of  any  discontinuity  employed  the 
frozen  flow  model.  The  frozen  flow  model  does  not  account 
for  chemical  reactions  within  the  nozzle.  This  frozen  flow 
assumption  was  justified  by  Snelling's  finite  rate  chemistry 
investigations,  which  showed  that  the  rapid  expansion  within 
the  hypersonic  nozzle  essentially  froze  the  chemical 
reactions  within  the  nozzle  (15:4-4).  However,  the  static 
temperatures  in  the  flows  considered  by  Snelling  were  only 
2000  K  for  uniform  flow  profile  studies,  and  the  nozzle 
inlet  profiles  considered  in  this  study  approached  3000  K. 
The  finite  rate  chemistry  effects  become  increasingly  more 
important  as  the  temperature  increases  (1:373-375). 
Therefore,  SCHNOZ  code  runs  were  also  performed  using  finite 
rate  chemistry  kinetics  in  the  flowfield.  Additionally, 
perfect  gas  runs  were  made  using  the  efficient  upwind  code 
to  demonstrate  the  ability  of  the  perfect  gas  model  to 
capture  the  trends  in  the  hypersonic  nozzle  flowfield. 

The  initial  internal  flowfield  of  the  nozzle  was 
generated  from  a  RAMJET  Performance  Analysis  cycle  analysis 
computer  code  (11).  The  internal  flowfield  cycle  analysis 
solution  required  the  input  of  the  freestream  properties  at 
each  flight  condition.  The  freestream  properties  used  were 
determined  for  a  typical  hypersonic  flight  trajectory  with  a 
1000  psf  dynamic  pressure  loading  (8) . 

The  analysis  of  the  nozzle  CFD  code  results  was  made  at 
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three  points  along  the  1000  psf  flight  trajectory.  Flight 
conditions  of  Mach  10,  Mach  15,  and  Mach  20  were  considered. 
The  desired  goal  was  to  determine  the  most  effective  and 
efficient  CFD  code  for  use  across  the  entire  flight 
trajectory. 

Significant  areas  of  interest  in  comparison  of  the  CFD 
code  results  were  the  required  CPU  time,  the  code  accuracy 
in  capturing  the  flow  physics,  and  the  variation  of  code 
results  for  different  flight  conditions  and  different  nozzle 
entrance  profiles.  Additionally,  since  the  SCHNOZ  code 
proved  to  be  very  sensitive  to  grid  spacing  (13:4-2)  the 
efficiency  and  accuracy  of  grid  clustering  techniques  were 
investigated.  The  AFIT  Computer  Lab  facilities  were  used  to 
conduct  the  computational  analysis  and  comparison  with  all 
computations  being  performed  using  double  precision  floating 
point  arithmetic  64-bit  SPARCstation  2  machines. 

1 . 4  Approach 

Running  the  SCHNOZ  code  in  the  inviscid  mode  allowed 
for  a  comparison  between  the  FDS  Godunov-based  scheme  and 
MacCormack's  scheme  in  the  solution  of  the  steady  Euler 
equations.  Prior  to  the  actual  nozzle  flow  analysis,  a 
determination  of  the  trends  exhibited  by  each  code  was  made 
by  comparing  the  results  of  the  FDS  Godunov  scheme  and  the 
modified  MacCormack  scheme  to  a  more  simple  flow  situation; 
the  oblique  shock  wave.  The  oblique  shock  wave  has  an  exact 
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analytical  solution  which  may  be  compared  to  the  results  of 
the  two  different  schemes.  Comparison  of  the  scheme  results 
to  the  exact  solution  indicated  the  ability  of  each  method 
to  capture  the  true  physics  of  the  flow. 

Because  of  the  inherent  non-linearity  of  the  governing 
equations  and  the  complicated  interactions  within  the 
nozzle,  no  known  analytical  solution  exists  for  the 
hypersonic  nozzle.  It  was  expected  that  the  full  viscous, 
turbulent,  PNS  equations  modeled  in  the  SCHNOZ  program  would 
provide  more  accurate  information  on  the  propulsive  nozzle 
flowfield  as  it  accounted  for  viscous  interactions  as  well 
as  a  turbulence  model.  For  this  reason  the  viscous  SCHNOZ 
code  was  used  as  the  basis  for  accuracy  comparisons  between 
the  codes.  Investigations  of  each  flight  condition  were 
made  operating  the  SCHNOZ  code  in  both  viscous  and  inviscid 
modes.  The  CPU  time  required  for  each  code  solution  was 
measured  and  used  to  determine  the  efficiency  of  the  two 
codes . 

Investigations  at  each  flight  trajectory  point  were 
made  for  an  isolated  nozzle  with  no  external  flow 
interactions.  This  allowed  for  a  clean  comparison  between 
the  two  codes  and  the  ability  of  each  to  resolve  the 
gradients  within  the  nozzle  flowfield.  For  a  combined 
(internal  with  external)  flowfield,  both  computer  programs 
used  non-physical  procedures  to  obtain  solutions  across  the 
contact  surface  downstream  of  the  engine  cowl.  The  FDS 
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program  obtained  a  solution  by  employing  two  solid  surface 
boundary  conditions  and  requiring  the  pressure  to  match 
across  the  surface  (6:193).  The  SCHNOZ  program  solution  is 
obtained  by  Implementing  an  embedded  grid  scheme  at  an 
abritrary  distance  downstream  of  the  cowl,  where  the  grid 
spacing  is  based  on  the  chemical  species  parameter  profile 
(4:909).  For  thrust  considerations,  the  isolated  nozzle 
flowfield  was  shown  by  Snelling  to  be  equivalent  to  the 
combined  flowfield  for  a  given  cowl  geometry  (15:4-10). 
Therefore,  the  performance  of  each  code  on  the  isolated 
nozzle  was  a  good  indicator  of  the  numerically  determined 
nozzle  thrust  for  each  configuration. 
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Figure  l.l.  Typical  hypersonic  vehicle  with  airframe- 
integrated  nozzle  (6:5) 
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II .  Theory 


2.1  Governing  Equations 

The  governing  equations  of  fluid  dynamics  and  their 
abilities  to  solve  fluid  flow  phenomena  are  summarized  in 
Figure  2.1  (9:659).  From  the  conditions  present  in  the 
hypersonic  propulsive  nozzle  flowfield  discussed  previously, 
the  Navier-Stokes  equations  are  necessary  to  fully  solve  the 
flowfield.  Full  Navier-Stokes  equations  incorporate  the 
viscous  effects  of  molecular  and  turbulent  dissipation  and 
diffusion  (9:659).  In  the  hypersonic  environment,  radial 
changes  in  the  viscous  terms  dominate  axial  or  streamwise 
changes  (1:344)  and  an  accurate  solution  to  the  flowfield 
may  be  obtained  by  dropping  the  streamwise  derivatives  of 
viscous  terms.  The  resulting  set  of  equations,  with  the 
elimination  of  all  time  derivatives,  are  termed  the 
parabolized  Navier-Stokes  (PNS)  equations.  Full  discussions 
of  the  Navier-Stokes  and  PNS  equations  and  their  numerical 
implementation  are  presented  in  Anderson,  Tannehill,  and 
Pletcher  (2) .  The  PNS  equations  employed  in  the  SCHNOZ 
computer  program  are  presented  in  Appendix  A. 

The  effects  of  viscosity  and  the  conduction  of  heat  are 
not  included  in  the  Euler  equations.  In  the  hypersonic 
nozzle  design,  however,  the  Euler  equations  are  a  valuable 
tool  as  the  determination  of  nozzle  thrust  is  dominated  by 
the  inviscid  pressure  effects.  With  increasing  velocities 
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and  the  favorable  pressure  gradient  in  the  expansion  of  the 
nozzle,  the  viscous  boundary  layer  growth  is  limited,  thus 
having  a  lessening  effect  on  the  overall  flow  phenomena 
(2:235-236).  The  increasing  velocity  resulting  from  the 
nozzle  expansion  also  acts  to  diminish  the  free  shear  layer 
interactions  within  the  nozzle  flow.  The  Euler  eguations  are 
capable  of  capturing  the  relevant  flow  phenomena  and 
discontinuities  within  the  hypersonic  propulsive  nozzle 
f lowf ield. 

2.2  Thermodynamic  Models 

The  thermodynamic  model  employed  greatly  influences  the 
flow  properties  downstream  of  discontinuities  within  the 
hypersonic  nozzle.  The  simplest  model  is  the  perfect  gas. 

A  perfect  gas  model  assumes  no  intermolecular  forces  and 
constant  specific  heats,  Cy  and  Cp,  resulting  in  a  constant 
specific  heat  ratio,  -y.  The  perfect  gas  obeys  the  thermal 
equation  of  state  (1:381): 

p  =  pRT  (2.1) 

The  perfect  gas  assumption  breaks  down  when  high 
temperatures  are  encountered  within  the  nozzle  f lowf ield. 

In  real  gases,  y  is  not  a  constant,  but  is  dependent  on  the 
temperature  of  the  flowfield.  A  calorically  imperfect,  but 
thermally  perfect  gas  model  allows  for  the  temperature 
dependence  of  the  flow  by  considering  the  internal  structure 
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of  the  gas.  The  calorically  imperfect  gas  assumes  that  the 
thermal  equation  of  state,  Eg.  (2.1),  still  holds,  but 
accounts  for  a  non-constant  7  resulting  from  increased  modes 
of  energy  available  to  the  gas  molecule;  see  Appendix  B. 

The  calorically  imperfect  gas  assumption  is  coupled  with  a 
frozen  flow  model  which  assumes  that  the  molecular 
composition  of  the  gas  remains  unchanged  throughout  the 
flowfield.  The  complex  chemical  mixtures  which  exit  the 
SCRAM jet  combustor  can  be  treated  as  a  mixture  of  thermally 
perfect  gases.  The  resulting  flowfield  properties  such  as 
enthalpy  and  specific  internal  energy  are  no  longer  simple 
functions  of  the  temperature,  but  rather  are  dependent  on 
the  chemical  composition  of  the  mixture  and  the  functional 
dependence  on  temperature  of  each  species  contained  within 
the  mixture. 

In  real  gases,  the  flow  properties  are  a  function  of 
not  only  temperature,  but  pressure  and  time  as  well  (7:99). 
The  nonequilibrium  chemically  reacting  gas  thermodynamic 
model  accounts  for  the  chemical  reactions  within  the 
flowfield.  The  nonequilibrium  chemically  reacting  gas 
allows  for  a  changing  molecular  composition  based  on  finite- 
rate  chemistry  kinetics.  As  discussed  in  Appendix  B,  a 
finite-rate  chemistry  kinetics  package  is  coupled  with  the 
governing  equations  of  the  flow  and  determines  the  chemical 
composition  and  resulting  flow  properties  based  on  the 
reaction  rates  for  a  given  set  of  species  and  chemical 


2-3 


reactions.  The  nonequi librium  chemically  reacting  gas  is  a 
more  complete  model  of  the  flow  physics,  but  also  more 
complex  to  incorporate  in  a  solution  procedure. 

2.3  Computational  Methodologies 

2 .3 .1  PNS  Solvers 

PNS  equation  solvers  are  in  widespread  use  in  the 
analysis  of  the  hypersonic  flow  regime  (1:360).  The 
popularity  of  the  PNS  solvers  is  attributable  to  the 
efficiency  in  predicting  the  hypersonic  flowfield  with  a 
great  savings  in  computational  storage  and  time  in 
comparison  to  a  full  Navier-Stokes  solution  (2:420). 
Unfortunately,  to  obtain  an  accurate  solution  requires  a 
greater  amount  of  user  interaction  and  adjustment  of  various 
input  parameters  (1:348-349).  The  parabolization  of  the 
Navier-Stokes  equations  allows  for  an  explicit  downstream 
marching  finite-difference  solution  technique. 

A  Scramjet  Hypersonic  Nozzle  (SCHNOZ)  parabolized 
Navier-Stokes  computer  code  was  developed  by  Science 
Applications  International  Corporation  (SAIC)  for  the  NASP 
Program  (20) .  The  SCHNOZ  Code  is  the  nozzle  analysis 
component  of  an  integrated  system  of  two-dimensional  PNS 
codes  for  analyzing  SCRAMjet  propulsive  nozzle  flowfields 
(20:7).  SCHNOZ  unifies  previous  work  in  rocket  propulsive 
nozzles  and  the  extension  of  the  PNS  equations  to  supersonic 
mixing  problems  and  finite  rate  chemistry  (20:4). 
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SCHNOZ  uses  the  MacCormack  method  to  explicitly  solve 
the  PNS  equations  in  conservative  form  (20:4).  The 
MacCormack  predictor-corrector  method  determines  the 
downstream  flux  values  based  on  the  average  of  a  predicted 
and  then  corrected  downstream  flux  derivative  (2:482-485). 

To  dampen  oscillations  inherent  in  the  MacCormack  scheme  in 
the  presence  of  strong  discontinuities,  a  self-adjusting 
hybrid  scheme  developed  by  Harten  and  Zwas  is  added  to 
modify  the  MacCormack  scheme  (20:24). 

To  account  for  chemically  reacting  flow  in  the  nozzle, 
SCHNOZ  employs  a  finite-rate  chemical  kinetics  package 
developed  by  AeroChem  (20:20-21).  Additionally,  the  SCHNOZ 
code  contains  two  high  Reynolds  Number  turbulence  models  to 
handle  the  turbulent  nature  of  the  nozzle  inlet  flow  exiting 
from  the  combustor  (20:9). 

The  PNS  equations  are  hyperbolic  in  nature,  thus 
requiring  supersonic  flow  for  a  valid,  well-posed  problem 
(2:24).  The  addition  of  the  viscosity  effects  causes  a 
small  subsonic  region  at  the  viscous  boundary  layer 
interaction  region  (1:204).  This  subsonic  region  is 
mathematically  elliptic  and  the  downstream  marching 
technique  of  MacCormack' s  method  is  invalid  in  this  region 
(1:204).  The  SCHNOZ  code  removes  this  subsonic  portion  of 
the  boundary  layer  and  enforces  a  slip  boundary  condition, 
resolving  the  flow  at  the  wall  based  on  a  viscous- 
characteristic  formulation  related  to  the  pressure  and  flow 
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angle  at  the  boundary  (20:13-14). 

Snelling  (15)  used  the  SCHNOZ  code  to  examine  the 
effects  of  non-uniform  nozzle  entrance  profiles  on  pitching 
moments.  Snelling  investigated  several  inlet  nozzle  profiles 
on  an  isolated  and  combined  internal  flow/external  flow 
nozzle.  In  that  study  Snelling  showed  that  the  chemical 
kinetics  considerations  are  unnecessary  extra  computational 
time,  producing  nearly  identical  results  to  the  frozen  flow 
model  (15:5-1).  Snelling  also  found  that  the  addition  of 
the  external  flow  interaction  had  little  effect  on  the 
overall  nozzle  thrust  (15:4-11). 

2.3.2  Flux’-Difference-Splitting 

The  solution  of  discontinuities  in  the  nozzle  flowfield 
is  representative  of  the  Riemann  problem.  The  Riemann 
problem  describes  the  collapse  of  flowfield  discontinuities 
to  a  local  point  of  consideration  (6:2,11-17).  These 
discontinuities  give  rise  to  fluxes  that  are  not  present  in 
the  initial  value  line.  The  differences  between  these 
fluxes  and  the  initial  values  of  the  flowfield  are  split 
along  the  characteristic  waves  developed  by  the  local 
collapse  of  the  discontinuity  (6:2).  This  technique  is 
known  as  the  flux-difference-splitting  (FDS)  method.  Full 
details  on  the  FDS  method  and  the  solution  of  the  Riemann 
problem  using  a  Godunov  scheme  are  given  in  reference  (6) 
and  summarized  in  Appendix  C. 
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The  MacCormack  scheme  applied  to  the  Riemann  problem 
and  FDS  was  analyzed  by  Steger  and  Warming  (17)  for  the 
solution  of  time  dependent  flow  in  a  shock  tube.  Similarly, 
Roe  (12)  used  a  density  ratio  across  the  discontinuities  to 
influence  the  fluxes  considered  in  the  Riemann  problem  to 
solve  the  unsteady  Euler  equations.  Grossman  (7)  modified 
the  upwind  FDS  schemes  developed  by  Van  Leer,  Steger  and 
Warming,  and  Roe  to  incorporate  a  real-gas  equivalent  ratio 
of  specific  heats  based  on  equilibrium  chemical  reactions  in 
high  temperature  flow.  Sod  (17)  provides  an  excellent 
survey  of  the  different  solution  schemes  applied  to  the 
Riemann  problem  for  the  unsteady  Euler  equations. 
Chakravarthy  (3)  provides  a  comprehensive  comparison  of 
upwind  schemes  based  on  the  FDS  method  for  unsteady 
flowfields  that  are  integrated  in  time  to  the  steady  state. 

However,  the  determination  of  the  nozzle  thrust  for  a 
propulsive  vehicle  requires  only  the  steady-state  solution 
of  the  flowfield  (6:2).  The  use  of  the  time-dependent 
schemes  mentioned  above  are  not  very  efficient  when  advanced 
to  the  steady  state.  It  is  more  efficient  to  directly  solve 
the  steady  form  of  the  Euler  equations.  Pandolfi  (10), 
using  a  Godunov  scheme,  applied  the  unsteady  FDS  method  to 
the  solution  of  steady  supersonic  flows. 

Doty  (6)  employed  an  upwind  FDS  method  based  on  the 
Godunov  scheme  to  solve  the  Riemann  problem.  Doty's  work 
showed  that  the  shock  capturing  ability  of  the  first-order 
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accurate  upwind  FDS  method  provided  solutions  comparable  in 
accuracy  to  a  higher-order  MacCormack  method,  but  with  a 
tremendous  savings  in  computational  time  (6:75) . 

Herring  (8)  demonstrated  the  utility  of  the  upwind  FDS 
code  developed  by  Doty  in  an  optimization  study  of  nozzle 
contours  for  a  typical  hypersonic  flight  trajectory  assuming 
a  calorically  perfect  gas  in  the  nozzle.  Schieve  (13) 
modified  Doty's  code  to  use  a  calorically  imperfect 
(thermally  perfect)  gas  model  and  showed  a  difference  in 
nozzle  thrust  of  over  16%  compared  to  the  perfect  gas  model. 
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Figure  2.1.  Capabilites  of  gas  dynamics  equations  (9:659). 
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III.  Methodology 


3.1  FDS  Code  Description 

3.1.1  Governing  Flow  Equations 
The  FDS  code  of  Doty  (6)  solves  the  following 
divergence  vector  form  of  the  Euler  equations: 


dE 

dx 


+ 


dF 

dy 


=  0 


where  the  E  and  P  flux  vectors  are  given  by 


(3.1) 


pu 

pv 

pu^+p 

F  = 

pvu 

puv 

0  * 

pv2+p 

u{pe+p). 

v(pe+p). 

The  vector  components  represent  the  continuity,  x  momentum, 
y  momentum,  and  energy  equations,  respectively.  The 
governing  vector  equation  is  independent  of  the 
thermodynamic  models  considered  in  the  FDS  study  and  is 
applicable  to  both  perfect  and  calorically  imperfect  gas 
flows. 


3.1.2  Computational  Scheme 

The  governing  equations  are  transformed  from  physical 
space  to  the  computational  space  assuming  the  following 
transformation : 


The  governing  equation. 


5 

T|  =  T)  (x,y) 

Eq.  (3.1),  is  transformed  in 


(3.3) 
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computational  space  as: 


d(E)  _  ^  dim  _  „  dim 
~dr  9ri  dti 


(3.4) 


Detailed  information  on  the  coordinate  transformation  is 
presented  in  reference  (6:142-143). 

The  transformed  governing  equations  are  numerically 
solved  in  computational  space  using  a  first-order  accurate 
FDS  method.  The  FDS  method  developed  by  Doty  is  based  on 
the  Godunov  initial-value  Riemann  problem.  The  Riemann 
problem  and  the  first-order  accurate  upwind  FDS  method  are 
described  in  detail  in  Appendix  C.  Essentially,  the  Riemann 
problem  is  used  to  resolve  discontinuities  in  the  flowfield 
into  numerical  fluxes.  The  resulting  Riemann  fluxes  are 
propagated  along  characteristics  determined  by  the  local 
wave  angle  (6:9) . 

To  axially  march  the  solution  from  station  i  to  i+1  in 
Figure  3.1,  the  Godunov-based  FDS  method  resolves  the 
Riemann  fluxes  and  then  computes  flux  differences  based  on 
the  initial  value  flux  at  axial  location  i.  These  flux 
differences  are  then  split  based  on  their  direction  of 
propagation  such  that  the  resulting  flux  difference  is  sent 
in  the  physically  correct  direction. 

The  stencil  for  the  first-order  accurate  FDS  method  is 
presented  in  Figure  3.1.  Note  in  Figure  3 . 1  that  the 
Riemann  fluxes  are  resolved  at  the  grid  midpoints,  j+1/2  and 
j-1/2.  The  first-order  accurate  upwind  FDS  scheme  uses 
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positively  biased  information  from  node  j-1/2  and  negatively 
biased  information  from  node  j+1/2  to  advance  the  solution 
from  node  (i,j)  to  {i+l,j)  (13:17).  The  first-order 
accurate  FDS  finite-difference  approximation  to  the 
transformed  governing  equation.  Eg.  (3.4),  for  an  interior 
grid  point  is  given  by  (6:15): 

-  A5Ti,{dE5.i/2  +  dEj.^/^]  -  A§Ti3,{dP5.i/2  +  (3.5) 

It  should  be  noted  that  the  FDS  finite-difference  equation. 
Eg.  (3.5),  is  independent  of  the  Riemann  problem  solution 
procedure  and  the  thermodynamic  model. 

3.1.3  Boundary  Conditions 

The  interior  grid  point  solution,  Eg.  (3.5),  is  not 
valid  at  boundaries,  as  it  would  require  information  outside 
the  physical  domain.  As  shown  in  Figure  3.2,  there  can  be 
no  negatively  biased  information  from  a  node  j+1/2 
influencing  the  solution  because  of  the  solid  wall  upper 
boundary.  The  solution  procedure  for  the  FDS  solid  wall 
boundary  is  a  two  step  wave-corrector  requiring  the  velocity 
vector  to  be  tangent  to  the  solid  boundary  (6:190).  A 
contact  surface  boundary  point  ,  such  as  would  exist  for  the 
interaction  between  the  internal  nozzle  flow  and  external 
freestream,  essentially  uses  a  coupled  solid  wall  boundary 
condition.  An  iterative  solution  procedure  is  necessary  to 
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match  the  pressure  and  flow  angle  for  the  contact  contact 
surface  (6:193).  Details  on  the  first-order  accurate  FDS 
boundary  point  calculations  are  presented  in  reference 
(6:189-194) . 


3.1.4  Decoding  the  Solution 

Once  the  solution  has  been  marched  downstream  to  the 
next  axial  location,  i+1,  the  primitive  flow  variables  must 
be  extracted  from  the  newly  determined  E  and  P  flux  vectors 
so  that  the  Riemann  problem  may  be  solved.  This  allows  the 
solution  to  continually  march  downstream  in  the  axial 
direction.  The  extraction  of  the  primitive  varaibles  from 
the  E  and  P  flux  vector  solutions  is  referred  to  as  decoding 
the  solution.  The  decode  procedure  differs  depending  on  the 
thermodynamic  model  employed  in  the  calculation.  The 
perfect  gas  decode  procedure  is  a  closed  form  solution 
procedure  based  on  the  flux  vector  components  and  the 
conservation  of  stagnation  enthalpy.  The  calorically 
imperfect  gas  decode  procedure  is  an  iterative  process  based 
on  temperature  that  must  incorporate  the  imperfect  gas  model 
into  the  conservation  of  stagnation  enthalpy.  Differences 
in  the  decode  procedure  for  the  two  different  thermodynamic 
models  are  presented  in  Appendix  D. 
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3.2  PNS  Code  Description 

3.2.1  Governing  Flow  Equations 

The  SCHNOZ  code  solves  the  following  conservative  form 
of  the  PNS  equations: 


dE  ^  ^ 
dx  dy 


(3.6) 


where  the  flux  vectors  are 
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The  vector  components  of  Eq.  (3.7)  represent  the  continuity, 
X  momentum,  y  momentum,  energy  equation,  and  species 
continuity  equations,  respectively.  The  governing  vector 
equation  is  not  independent  of  the  thermodynamic  models 
considered  in  the  SCHNOZ  study.  The  species  continuity 
equation  is  necessary  for  finite-rate  chemistry 
considerations.  For  perfect  or  calorically  imperfect  gas 
frozen  flow  studies,  the  species  continuity  equation  is  not 
implemented.  Additionally,  for  inviscid  flow  calculations 
of  a  perfect  or  calorically  imperfect  frozen  gas,  the 
viscous  terms  are  dropped  and  the  PNS  equations  reduce  to 
the  Euler  equations. 

The  governing  equations  are  transformed  from  physical 
space  into  rectangular  computational  coordinates,  |  and  tj. 
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through  the  following  transformation: 

n  =  [y  -  /  [Vaix)  -  y^(x)] 

where  yu  and  Yl  represent  the  upper  and  lower  physical 
boundaries  respectively  (20:21). 

3.2.2  Computational  Scheme 

The  transformed  governing  equations  are  numerically 
solved  in  the  computational  domain  using  the  explicit  two- 
step  MacCormack  algorithm.  The  MacCormack  algorithm 
spatially  marches  the  solution  from  an  initial  data  line,  i, 
downstream  to  a  new  axial  location,  i+1.  The  first  step  of 
the  MacCormack  method  calculates  a  downstream  predicted 
value  of  the  E  flux  vector,  denoted  by  E*  ,  using  the 
upstream  values  of  the  E,  P,  6,  and  f  flux  quantities 
(1:197).  The  second  step  of  the  MacCormack  method  then 
uses  the  predicted  flux  values,  B*,  P*,  6*,  and  f’,  to  correct 
the  downstream  flux  quantities  (1:198).  The  SCHNOZ  program 
employs  the  following  second-order-accurate,  central- 
difference  predictor  and  corrector  formulas  at  an  interior 
grid  point,  (i,j),  to  advance  the  solution  from  an  axial 
location  i  to  i+1  (15:3-2): 
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Predictor  Step: 


=  Bj  -  -  {l-2e)Fj  -  eFj.j^  - 

,  ^  -  f,  -  A-  {£,  -  f,.,) } 


(3.9) 


Corrector  Step: 


=  l[E^j  +  Ej  -  +  (l-2e)Fj  +  {e-DF^j.^} 

-a*jA5  +  [A**  -  f*j)  -  A*  {£*J  -  :^-.i)}] 

Ax\^ 


(3.10) 


where 


A*  = 


and  for  f  =  u,  v,  H,  and  aj 


o 


U 


0^=1 


= 


=  PR 


(3.11) 


(3.12) 


The  e  term  in  Eqs.  (3.9)  and  (3.10)  is  an  alternating  switch 
(e  =  0  at  even  steps;  e  =1  at  odd  steps)  used  to  provide 
nonbiased  convective  differencing  (20:23). 

The  basic  MacCormack  scheme  has  been  shown  to  overshoot 
and  undershoot  flow  properties  downstream  of  flowfield 
discontinuities  (20:23,  2:146-147,  and  6:24,31).  The 
oscillation  of  flow  properties  about  a  discontinuity  is 
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remedied  using  the  self  adjusting  hybrid  scheme  of  Harten 
and  Zwas  (20:23).  The  Harten  and  Zwas  hybrid  scheme  adds  a 
dissipative  term  to  the  conservative  flow  equations  of  the 
MacCormack  scheme  based  on  the  pressure  at  surrounding  grid 
points  (20:23).  From  a  corrected  value  of  the  flux  variable 
E,  the  addition  of  a  viscous  term  to  the  MacCormack 
solution  is  implemented  to  yield  the  self-adjusting  hybrid 
value,  E,  using  the  following  finite-difference 
representation  (20:23): 


(3.13) 


The  d  terms  are  dimensionless  parameters  varying  from  0  to 
1,  based  upon  the  pressure  at  the  surrounding  grid  points  to 
provide  appropriate  levels  of  damping  to  the  solution 
(20:23) . 


3.2.3  Boundary  Conditions 

The  SCHNOZ  program  uses  the  forward  or  backward 
MacCormack  algorithm  with  an  Abbett  correction  procedure  to 
handle  solid  wall  boundary  conditions  (20:31).  The  Abbett 
correction  procedure  enforces  the  surface  flow  tangency 
condition  at  the  wall.  A  corrected  pressure,  Pj,  is 
determined  at  the  wall  by  considering  the  wall  to  be 
composed  of  a  series  of  infinitesimal  simple  expansion  or 
compression  waves  for  a  small  change  in  surface  slope,  A6 . 
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The  corrected  pressure  is  found  from  (15:3-3); 


Pc  =  p  -  (3.14) 

where  p  and  M  in  Eq.  (3.14)  are  the  values  at  the  boundary 
determined  using  the  MacCormack  algorithm.  The  corrected 
pressure  is  then  used  to  correct  the  density  and  the  axial 
and  radial  velocity  components  are  determined  from  the 
actual  slope  of  the  wall.  More  details  are  given  by  Wolf  et 
al.  (20:33-34). 

3.2.4  Turbulence  Model 

The  SCHNOZ  program  contains  two  variants  of 
turbulence  modelling  equations,  the  ke  and  kW  high  Reynolds 
number  versions.  The  turbulence  model  used  in  this  study  is 
the  common  two-equation  k€  eddy  viscosity  model  which  models 
the  turbulent  kinetic  energy,  k,  and  energy  dissipation,  e, 
detailed  in  Appendix  A.  The  ke  model  was  chosen  because  the 
kW  model  has  had  difficulties  in  applications  to  wall 
bounded  shear  flows  such  as  those  in  a  nozzle  (5:507). 
Additionally,  the  implementation  of  the  ke  model  in  the 
SCHNOZ  code  is  determined  by  a  correlation  to  experimental 
jet  mixing  results  (5:507). 

3.2.5  Decoding  the  Solution 

Unlike  the  FDS  algorithm,  the  primitive  variables  are 
not  a  computational  necessity  at  every  axial  location  in  the 
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downstream  marching  of  the  SCHNOZ  program.  Advances  from 
one  axial  position  to  the  next  using  the  MacCormack  method 
are  based  only  on  the  flux  vectors.  The  extraction  of 
primitive  variables  from  the  flux  vectors  (decoding  the 
solution)  is  nec  ssary  however,  to  determine  the  artificial 
damping  of  the  Harten  and  Zwas  hybrid  scheme,  which  is 
dependent  on  the  pressure  at  surrounding  nodes.  The  decode 
procedure  is  also  used  to  obtain  the  primitive  variables  at 
desired  axial  locations.  The  SCHNOZ  program  uses  a 
different  decode  procedure  based  on  the  thermodynamic  model 
employed.  For  the  perfect  gas  model,  a  simple  closed  form 
decode  procedure  is  employed  based  on  the  flux  vector 
components.  For  the  thermally  perfect  frozen  flow  and  the 
finite-rate  chemistry  thermodynamic  models,  a  different  real 
gas  decode  procedure  is  employed  depending  on  whether  the 
flow  is  supersonic,  M>1,  or  hypersonic,  M>5. 

For  supersonic  flows,  an  iterative  procedure  is 
employed  that  uses  an  assumed  pressure  value,  p' ,  to 
calculate  the  flow  velocity.  The  total  enthalpy  and 
velocity  are  then  used  to  calculate  the  enthalpy,  h' : 

h'  =  H'  -  ^  (u^  +  v^)  (3.15) 

Calculating  the  species  mole  fractions,  a,,  and  the 
temperature,  T',  from  the  assumed  value  of  p',  the  enthalpy 
is  also  calculated  using  the  species  formulation: 


3-10 


h*  =  Za^h^iT') 


(3.16) 


The  assumed  pressure,  p' ,  is  iterated  upon  until  the 
calculated  enthalpy  values  of  h'  and  h*  match  to  a 
prescribed  tolerance  (20:30). 

The  iteration  procedure  described  above  breaks  down  for 
the  expanding  hypersonic  nozzle  flow  where  the  kinetic  flow 
energy,  H,  is  orders  of  magnitude  greater  than  the  thermal 
contribution,  h  (20:31).  SCHNOZ  uses  an  iterative  decode 
procedure  for  hypersonic  flow  conditions  that  is  based  on  a 
linear  relation  between  the  axial  velocity  component  and 
temperature  (20:75) . 

3.3  Oblique  Shock  Validation 

3.3.1  Perfect  Gas 

A  demanding  test  of  the  two  numerical  methods  is  the 
oblique  shock  reflection  which  contains  large 
discontinuities  in  the  flowfield  properties  downstream  of 
the  shock  wave.  The  geometry  of  this  flow  situation  is 
shown  in  Figure  3.3.  Doty,  using  a  perfect  gas  assumption, 
compared  the  second-order  centered  MacCormack  algorithm  to 
the  first-order  upwind  Godunov  based  FDS  algorithm  using  a 
10  degree  ramp  shock  wave  geometry  with  a  Mach  2.2  incoming 
flow  (6:31).  Figure  3.4  shows  the  static  pressure 
distribution  along  the  top  wall  obtained  with  the  two 
different  methods  compared  to  the  exact  analytical  solution. 
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Doty  found  that  the  second-order  accurate  MacCormack  method 
did  no  better  than  the  first-order  accurate  Godunov  based 
FDS  method  in  capturing  the  location  of  the  shock  (6:24). 
Additionally,  the  MacCormack  method  overshot  the  exact 
solution  before  reaching  the  correct  pressure  rise,  while 
the  FDS  method  monotonically  approached  the  correct  pressure 
(6:24). 

The  MacCormack  method  implemented  by  Doty  did  not 
incorporate  the  Marten  and  Zwas  hybrid  scheme  to  remedy  the 
oscillations  downstream  of  the  discontinuity.  Sod  showed 
that  the  addition  of  the  artificial  viscosity  of  the  Marten 
and  Zwas  hybrid  scheme  smoothed  out  the  MacCormack  method  in 
the  area  of  discontinuities  while  retaining  second-order 
accuracy  in  the  smooth  portions  of  the  flow  (16:24).  The 
oblique  shock  reflection  geometry  verified  the  ability  of 
the  first-order  Godunov-based  FDS  scheme  to  capture  the 
physics  of  a  strong  flowfield  discontinuity.  This  provided 
the  initiative  to  implement  the  FDS  method  in  the  analysis 
of  the  hypersonic  propulsive  nozzle  flowfield. 

3.3.2  Imperfect  Gas 

The  MacCormack  method  as  implemented  in  the  SCMNOZ 
computer  program  incorporates  the  Marten  and  Zwas  hybrid 
damping  scheme.  A  comparison  between  the  imperfect  gas  FDS 
method  and  the  SCMNOZ  frozen  flow  inviscid  solution  was  made 
on  the  10  deg  oblique  shock  geometry  using  the  same  Mach  2.2 
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inlet  conditions  as  Doty  used  to  compare  the  perfect  gas  FDS 
and  MacCormack  methods.  To  allow  for  a  comparison  to  an 
exact  analytical  solution,  the  imperfect  gas  was  assumed  to 
be  Argon-free  air  in  the  standard  mole  fraction  ratios  of 
0.79  Nj  and  0.21  Oj.  The  analytical  solution  was  determined 
using  an  imperfect  gas  oblique  shock  wave  solver  provided  by 
Schieve  (13:91-95). 

Figure  3 . 5  shows  the  static  pressure  distribution  along 
the  upper  wall  obtained  with  the  two  different  methods 
compared  to  the  exact  analytical  solution.  The  FDS 
imperfect  gas  method  and  the  SCHNOZ  frozen  flow  solution 
equally  capture  the  the  location  of  the  shock  and  the  static 
pressure  downstream  of  the  flowfield  discontinuity. 
Additionally,  with  the  Marten  and  Zwas  damping  scheme  the 
second-order  MacCormack  method  no  longer  overshoots  the 
exact  solution  and  monotonically  approaches  the  correct 
pressure  as  does  the  first-order  accurate  Godunov-based 
scheme . 

3.4  Isolated  Hypersonic  Nozzle 

3.4.1  Freestream  Condi tions 

The  freestream  conditions  used  in  this  study  are  those  for  a 
typical  hypersonic  flight  trajectory  with  a  constant  dynamic 
pressure,  q„,  of  1000  psf.  Freestream  properties  were 
determined  at  Mach  numbers  of  10,  15,  and  20  along  this 
trajectory.  Assuming  the  freestream  conditions  obey  perfect 
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gas  relations,  the  freestream  static  pressure  is  found  from 
the  definition  of  Mach  number  and  the  constant  dynamic 
pressure : 


P  = 


2g, 

YACf 


(3.17) 


The  flight  trajectory  altitude,  z,  is  interpolated  from  the 
1962  Standard  Atmosphere  Table  assuming  p  =  p(z) .  The 
remaining  freestream  properties  are  calculated  using  the 
1962  atmosphere.  The  freestream  conditions  used  in  this 
study  are  tabulated  in  Table  3.1. 


3,4.2  Internal  Flow  Conditions 

The  initial  conditions  for  the  isolated  nozzle 
investigations  were  generated  from  a  RAMJET  Performance 
Analysis  (RJPA)  cycle  analysis  code  (10) .  The  RJPA  code 
simulated  a  Supersonic  Combustion  RAMJET  (SCRAM jet)  engine 
operating  at  each  flight  condition  along  the  hypersonic 
trajectory.  The  SCRAM jet  was  assumed  to  operate  in  a  shock- 
on-lip  condition  with  the  maximum  capture  of  airflow  without 
spillage  (19:3).  The  shock-on-lip  condition  is  illustrated 
in  Figure  3 . 6  where  the  bow  shock  wave  emanating  from  the 
forebody  compression  rests  on  the  engine  inlet. 

The  capture  area  of  the  SCRAM jet  inlet  is  driven  by 
airframe  geometries  and  the  resulting  shock  structure 
upstream  of  the  inlet  (19:3).  For  this  study,  an  effective 
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bow  shock  of  8  degrees  was  used  to  simulate  a  2  degree  angle 
of  attack  and  6  degree  forebody  compression  on  the  underside 
of  the  hypersonic  vehicle.  The  SCRAM jet  inlet  and  capture 
areas  are  then  determined  at  each  flight  condition  for  the  8 
degree  compression. 

The  RJPA  program  was  run  simulating  an  air-breathing, 
hydrogen-fueled,  SCRAM jet  engine  cycle  using  a  14  species 
gas  mixture  with  some  frozen  Hj  carried  through  a  constant 
area  combustion  process  to  simulate  combustion  inefficiency. 
The  RJPA  program  calculates  the  following  flow  properties  at 
several  stations  in  the  SCRAMjet  engine:  T,  p,  p,  V,  7^,  Sf, 
Mf;  the  f  subscript  denotes  frozen  flow  values.  The  RJPA 
program  also  calculates  the  molecular  weight  (MW)  of  the 
mixture  and  the  mole  and  mass  fractions  of  each  species  in 
the  mixture  at  specified  engine  stations.  The  combustor 
exit  properties  calculated  by  the  RJPA  program  provide  the 
initial  conditions  for  the  internal  nozzle  studies. 

The  FDS  code  requires  an  input  value  of  the  gas 
constant,  R,  which  was  determined  using  the  universal  gas 
constant,  Rq,  and  the  MW. 


This  value  of  R  is  used  as  the  initial  condition  for  both 
perfect  and  imperfect  gas  frozen  flow  analyses  and  finite- 
rate  chemistry  flows.  The  MW  and  R  values  remain  constant 
for  the  perfect  gas  and  frozen  flow  models  as  the  molecular 
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composition  remains  fixed.  The  finite-rate  chemistry  flow 
allows  for  chemical  reactions  and  hence  a  change  in  MW  and 
the  gas  constant,  R.  The  ratio  of  specific  heats,  jf, 
determined  from  the  RJPA  analysis  is  used  as  an  initial 
condition  of  the  internal  nozzle  for  all  flows  considered. 
The  differences  between  the  perfect  gas,  calorically 
imperfect  frozen  gas,  and  finite-rate  chemistry  assumptions 
do  not  effect  the  initial  conditions,  but  become  apparent  in 
the  downstream  expansion  of  the  flow  in  the  nozzle;  see 
Appendix  B.  The  nozzle  inlet  parameters  for  the  FDS  code 
are  tabulated  in  Table  3.2  for  each  flight  condition  on  the 
hypersonic  flight  trajectory. 

The  internal  nozzle  inlet  conditions  for  the  SCHNOZ 
code  are  presented  in  Table  3.3  for  the  same  flight 
trajectory  points.  The  SCHNOZ  code  using  the  perfect  gas 
assumption  requires  not  only  R,  but  the  MW  of  the  mixture. 

In  essence,  the  mixture  of  gases  in  the  nozzle  is  treated  as 
a  perfect  gas  with  a  MW  equivalent  to  the  MW  of  the  14 
species  gas  mixture  calculated  by  the  RJPA  program. 

3.4.3  Nozzle  Configuration 

As  shown  in  Figure  3.7  (6:46),  the  hypersonic  nozzle 
wall  consists  of  two  sections.  The  first  section,  A-B,  is  a 
circular  arc  of  radius  r,  followed  by  a  parabolic  section, 
B-C,  with  an  attachment  angle,  9^.  For  the  isolated  nozzles 
used  in  this  study  the  lower  boundary,  or  cowl  wall,  E-F,  is 
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a  straight  wall  section. 

The  parabolic  nozzle  contour  used  for  this  study  is  the 
optimal  nozzle  contour  for  the  entire  hypersonic  flight 
trajectory  determined  from  the  FDS  perfect  gas  analysis 
performed  by  Herring  (9) .  This  analysis  arrived  at  a 
parabolic  nozzle  contour  with  an  attachment  angle  of  20.6 
degrees.  Additionally,  a  much  higher  circular  arc  attachment 
angle,  6^  =  38  deg,  is  considered  to  examine  the  ability  of 
each  code  to  compute  a  flowfield  with  a  shock  wave  (6:49). 
The  nozzle  geometries,  x/h  and  y/h,  are  characteristic  of 
geometries  considered  for  hypersonic  propulsive  nozzles 
(6:3).  The  nozzle  geometry  parameters  for  the  two  nozzles 
are  listed  in  Table  3.4. 

The  nozzle  geometry  utilized  in  the  SCHNOZ  code  is 
flipped  upside  relative  to  the  FDS  nozzle  geometry,  however 
the  pertinent  specifications  remain  the  same.  The  x  and  y 
coordinates  for  each  wall  location  are  determined  from  the 
FDS  analysis  parabolic  contour  generator  and  are  used  to 
generate  the  nozzle  wall  input  to  the  SCHNOZ  code.  The 
SCHNOZ  wall  contour  is  then  calculated  from  the  x  and  y 
coordinates  using  a  cubic  spline  interpolation  routine.  The 
cubic  spline  interpolation  can  lead  to  small  differences 
between  the  two  codes  in  the  x  and  y  coordinates  and  nozzle 
wall  angles  at  each  wall  location. 
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Table  3.1.  Freestream  conditions  for  typical  hypersonic 
flight  trajectory 


Free  stream  parameter 

i 

1 

0 

CM 

It 

8 

s 

altitude,  (km) 

33.791 

39.581 

43.934 

static  pressure,  p  (N/m^) 

684.0 

304.0 

171.0 

static  temperature,  T  (K) 

233.2 

249.2 

261.2 

density,  p  (kg/m^) 

0.01022 

0.00425 

0.00228 

velocity,  V  (m/s) 

3061.4 

4746.9 

6479.9 

specific  heat  ratio,  7 

1.4 

1.4 

1.4 

gas  constant,  R  (J/kg/K) 

287.06 

287.06 

287.06 

Table  3.2.  Nozzle  inlet  conditions,  FDS  inputs 


Internal  Flow  Parameter 

=  10 

mmm 

S 

8 

II 

to 

0 

Mach  number 

2.0582 

3.4736 

4.6791 

static  pressure,  p  (N/m^) 

373,566 

199,752 

133,921 

static  temperature,  T  (K) 

2986.1 

3062.9 

3176.3 

specific  heat  ratio,  7 

1.25219 

1.25966 

1.27186 

gas  constant,  R  (J/kg/K) 

355.9 

365.7 

381.2 

Species  Mass  Fraction  for 

Frozen  Flow  | 

N2 

0.726652 

0.724947 

0.722695 

O2 

0.023287 

0.029532 

0.035884 

Ar 

0.012877 

0.012877 

0.012877 

H 

0.000701 

0.001445 

0.002840 

OH 

0.022336 

0.030736 

0.040664 

H2 

0.005085 

0.006512 

0.008066 

NO 

0.015743 

0.018386 

0.024192 

H7O 

0.188322 

0.164449 

0.132859 

0 

0.004994 

0.010106 

0.019907 
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Table  3.3.  Nozzle  inlet  conditions,  SCHNOZ  inputs 


Internal  Flow  Parameter 

M„=  10 

1 

1 

M„=  20 

Mach  number 

2.0582 

3.4736 

4.6791 

static  pressure,  p  (atm) 

3.6867 

1.9714 

1.3217 

static  temperature,  T  (K) 

2986.1 

3062.9 

3176.3 

axial  velocity,  u  (ft/s) 

7789.1 

13,537 

19,052 

radial  velocity,  v  (ft/s) 

0 

• 

0 

0.0 

0.0 

1  Gas  Mixture  Parameters  for  Perfect  Gas  Flow  | 

molecular  weight,  MW 

23.3635 

22.7327 

21.8060 

specific  heat  ratio,  7 

1.25219 

1.25966 

1.27186 

1  Species  Mole  Fraction  for  Frozen  and  Finite-Rate  Flow  | 

H 

0.016243 

0.032577 

0.061418 

0.058920 

0.073441 

0.087213 

HjO 

0.244184 

0.207456 

0.160752 

N2 

0.613426 

0.595442 

0.569337 

0 

0.007291 

0.014355 

0.027121 

OH 

0.030678 

0.041072 

0.052118 

O2 

0.017000 

0.014681 

0.024445 

NO 

0.012253 

0.020976 

0.017572 

Table  3.4.  Nozzle  geometry  specifications 


Circular  arc  attachment  angle,  (deg) 

20.6 

38.0 

Inlet  height,  h  (inches) 

1.0 

1.0 

Nozzle  circular  arc  radius,  r  (inches) 

1.0 

1.0 

Nozzle  wall  end,  x/h 

100.0 

100.0 

Nozzle  wall  end,  y/h 

±25.0 

±25.0 

Nozzle  wall  end  angle,  B^  (deg) 

±10.01 

0 

• 

CO 

+1 

Figure  3.1,  Stencil  for  first-order  accurate  upwind  FDS 
method  (6:195). 


Figure  3.2.  stencil  for  first-order  accurate  upper  solid 
wall  boundary  point  (6:195). 
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Figure  3.3.  Geometry  for  oblique  shock  wave  study  (6:30). 


Figure  3.4.  Static  pressure  distribution  along  the  top  wall 
for  shock  wave  reflection  study,  FDS  and 
MacCormack  methods  (6:31). 
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Figure  3.5.  static  pressure  distribution  along  the  top  wall 
for  shock  wave  reflection  study,  FDS  and 
SCHNOZ  frozen  flow  solutions. 


Figure  3.6.  Shock-on-lip  SCRAMjet  operation  (19:27). 
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IV .  Findings 


4 . 1  Case  Summary 

Comparisons  of  the  FDS  and  SCHNOZ  codes  were  made 
using  several  different  combinations  of  nozzle  geometries, 
flight  conditions,  and  thermodynamic  models.  Direct 
comparisons  of  the  FDS  and  SCHNOZ  codes  were  made  using  an 
inviscid  frozen  flow  assumption  for  two  different  parabolic 
nozzle  contours,  20.6  deg  and  38  deg,  at  three  different 
flight  conditions,  Mach  =  10,  15,  and  20,  along  a  typical 
hypersonic  trajectory.  An  investigation  of  the  SCHNOZ 
implementation  of  viscosity  terms  was  made  by  comparing 
inviscid  and  viscous  code  results  for  the  same  nozzle 
geometries  and  flight  conditions.  Comparisons  of  the  frozen 
flow  and  finite-rate  chemistry  models  in  the  SCHNOZ  code 
were  also  made  for  the  same  combinations  of  nozzle 
geometries  and  flight  conditions.  The  nozzle  geometries  and 
flight  conditions  were  also  used  to  compare  perfect  gas  and 
frozen  flow  FDS  results.  The  comparisons  are  presented  in 
plots  of  static  pressure  and  temperature  distributions  along 
the  nozzle  wall,  pressure  contour  plots,  and  tables  of  CPU 
time. 
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4.2  Grid  Refinement 

4.2.1  SCHNOZ 

The  SCHNOZ  code  proved  to  be  very  sensitive  to  grid 
fineness  in  resolving  the  changing  flowfield  properties 
along  the  nozzle  wall  (15:4-2).  With  a  coarse  grid,  51 
points  in  the  radial  direction,  the  MacCormack  algorithm  was 
unable  to  keep  track  of  the  rapidly  changing  flow  properties 
along  the  nozzle  wall  (15:4-2).  A  finer  grid,  101  radial 
points,  was  required  to  fully  resolve  the  complex  nozzle 
flowfield.  The  SCHNOZ  program  utilizes  equally  spaced 
radial  grid  intervals.  Ay.  The  allowable  axial  marching 
step.  Ax,  is  limited  by  a  combination  of  the  Courant- 
Friedrichs-Lewy  (CFL)  hyperbolic  stability  criterion  and  the 
parabolic  diffusion  criterion  for  viscous  flows  (20:32). 

4.2.2  FDS  Grid  Packing 

Grid  packing  techniques  are  available  in  the  FDS  code 
that  can  cluster  gridlines  toward  the  nozzle  wall  by 
reducing  the  radial  interval.  Ay.  The  grid  clustering  is 
dependent  on  the  packing  factor,  which  gives  the  fractional 
position  of  a  packed  grid  point  (6:147).  Figure  4.1  plots 
the  static  pressure  distribution  along  the  nozzle  wall  for 
the  FDS  perfect  gas  solution  using  a  grid  packing  value  of 
1.10  and  also  for  no  grid  packing.  Figure  4.2  plots  the 
static  pressure  distribution  along  the  nozzle  wall  for  the 
frozen  flow  FDS  solution  and  the  two  different  grid  packing 
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factors.  Figures  4.1  and  4.2  show  that  grid  clustering  is 
necessary  to  resolve  the  flowfield  properties  along  the 
nozzle  wall  without  oscillatory  behavior  using  the  FDS  code. 
Computations  were  made  with  101  radial  grid  points  for  both 
the  perfect  gas  and  frozen  flow  thermodynamic  models  on  the 
38  deg  parabolic  nozzle  contour  at  the  Mach  15  flight 
condition.  The  packing  factors  were  used  to  cluster  the 
gridlines  near  the  lower  and  upper  walls  to  capture  the 
rapidly  changing  properties  along  these  boundaries. 

Table  4 . 1  presents  a  comparison  of  the  number  of 
computational  planes  and  CPU  seconds  required  as  a  result  of 
different  grid  packing  factors.  The  CPU  seconds  are  for 
computations  using  double  precision  floating  point 
arithmetic  64-bit  SPARCstation  2  machines.  Because  of  the 
decrease  in  the  radial  interval,  the  grid  packing  causes  a 
corresponding  decrease  in  the  axial  step  size,  Ax,  to  meet 
the  CFL  stability  criterion.  No  significant  change  in  the 
wall  property  distributions  were  seen  with  the  increased 
grid  packing  for  a  packing  factor  of  1.05,  but  a  significant 
increase  in  CPU  time  was  required.  Attempts  to  use  51 
radial  grid  points  and  grid  packing  factors  were 
unsuccessful  as  the  FDS  code  failed  to  compute  a  non- 
oscillatory  solution  for  these  grids.  The  grid  packing 
factor  of  1.10  provided  a  non-oscillatory  solution  with  the 
least  inefficiency  in  terms  of  CPU  time  and  was  used 
throughout  the  remainder  of  the  study. 
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4.3  FDS  VS  SCHNOZ  Frozen  Flow  Comparisons 

Figures  4 . 3  through  4 . 14  present  a  comparison  between 
the  FDS  and  SCHNOZ  inviscid,  frozen  flow  solutions  to  the 
steady  Euler  equations.  The  comparison  between  the  two 
solutions  is  presented  in  contour  plots  of  static  pressure 
throughout  the  nozzle,  and  property  distributions  along  the 
nozzle  wall.  The  contour  plots  and  property  distribution 
plots  illustrate  the  particular  characteristics  of  the 
propulsive  nozzle  flowfield  and  demonstrate  the  ability  of 
each  computer  program  to  capture  the  flow  physics  in  the 
hypersonic  nozzle. 

The  static  pressure  contour  plot  for  the  FDS  imperfect 
gas  solution  with  the  20.6  deg  nozzle  attachment  angle  is 
illustrated  in  Figure  4.3  for  the  Mach  10  flight  condition. 
The  static  pressure  contour  plot  for  the  SCHNOZ  inviscid, 
frozen  flow  solution  for  the  same  flight  condition  and 
nozzle  geometry  is  illustrated  in  Figure  4.4.  Figures  4.3 
and  4.4  show  a  good  agreement  between  the  two  solution 
methods  for  the  general  characterization  of  the  nozzle 
flowfield.  Immediately  downstream  of  the  combustor  exit  the 
flow  undergoes  a  rapid  expansion  about  the  nozzle  attachment 
radius  (item  1) .  The  initial  expansion  is  reflected  off  the 
lower  boundary  and  continually  influences  the  flow 
downstream  in  the  nozzle.  As  the  rate  of  change  of  the 
slope  of  the  nozzle  wall  decreases,  a  recompression  is 
generated  in  the  flowfield  (item  2) .  This  recompression  is 
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reflected  off  the  lower  boundary  and  impinges  on  the  nozzle 
wall  (item  3) .  The  reflected  compression  is  then  carried 
throughout  the  nozzle  (item  4). 

The  calculated  pressure  distribution  along  the  nozzle 
wall  for  both  the  FDS  and  SCHNOZ  frozen  flow  solutions  is 
presented  in  Figure  4.5  for  the  Mach  10  flight  condition  and 
the  20.6  deg  nozzle  attachment  angle.  The  rapidly  expanding 
flow  around  the  nozzle  attachment  radius  causes  a  sharp 
initial  decrease  in  pressure.  The  pressure  rise  downstream 
at  an  X/h  of  between  5  and  10  results  from  the  recompression 
of  the  flow  due  to  the  change  in  wall  curvature  after  the 
parabolic  wall  attachment  point.  The  drop  in  pressure 
downstream  of  the  recompression  occurs  due  to  the  expanding 
nozzle  geometry  and  the  reflection  of  the  initial  expansion. 
The  increase  in  presssure  further  downstream  at  X/h=50  is  a 
result  of  the  reflected  recompression  wave  as  illustrated  in 
Figures  4.3  and  4.4.  The  expansion  and  pressure  drop 
gradually  proceeds  downstream  of  this  slight  recompression. 

The  temperature  distribution  along  the  nozzle  wall  for 
the  FDS  and  SCHNOZ  frozen  flow  solutions  for  the  Mach  10 
flight  condition  and  20.6  deg  nozzle  geometry  is  illustrated 
in  Figure  4.6.  The  temperature  distribution  along  the 
nozzle  wall  essentially  follows  that  of  the  pressure 
distribution.  The  temperature  decreases  in  regions  of 
expansion,  and  increases  across  the  recompressions. 

Figures  4.5  and  4.6  show  that  there  is  a  slight 
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difference  in  the  pressure  and  temperature  distributions 
along  the  nozzle  wall  for  the  two  different  code  solutions. 
The  recompression  pressure  peak  is  slightly  larger  in 
magnitude  and  located  slightly  further  downstream  for  the 
SCHNOZ  solutions.  The  pressure  and  temperature  rise  at  the 
reflected  recompression,  X/h=50,  are  more  pronounced  for  the 
SCHNOZ  solution  than  the  FDS  solution.  The  difference 
between  the  two  results  may  be  attributable  to  the 
additional  axial  steps  required  in  the  FDS  solution.  Table 
4.2  presents  a  comparison  of  the  number  of  computational 
planes  and  CPU  seconds  required  for  the  two  solutions.  For 
the  Mach  10  flight  condition  and  the  20.6  deg  nozzle 
attachmnet  angle,  the  FDS  solution  requires  almost  3  times 
as  many  computational  steps  (1461)  as  the  SCHNOZ  solution 
(562) .  The  increased  number  of  computational  steps  results 
from  the  grid  packing  utilized  in  the  FDS  solution  and  the 
CFL  stability  criterion  requirements.  The  increased  number 
of  computational  steps  adds  artificial  damping  to  the 
flowfield  solution,  thus  diminishing  the  pressure  and 
temperature  rise  at  the  compression  peaks.  Despite  these 
slight  differences  in  the  pressure  and  temperature  peaks, 
there  is  a  good  overall  agreement  between  the  two  code 
results. 

The  static  pressure  contour  plot  for  the  FDS  imperfect 
gas  solution  with  the  20.6  deg  nozzle  attachment  angle  is 
illustrated  in  Figure  4.7  for  the  Mach  15  flight  condition. 
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Figure  4.8  presents  the  same  information  for  the  SCHNOZ 
frozen  flow  solution.  Figures  4.7  and  4.8  show  a  good 
agreement  between  the  two  solution  methods.  The  rapid 
expansion  about  the  nozzle  attachment  circular  radius  occurs 
immediately  downstream  of  the  nozzle  inlet  (item  1) .  The 
change  in  nozzle  curvature  causes  a  recompression  (item  2) 
which  reflects  off  of  the  lower  boundary  (item  3)  and  is 
carried  downstream  (item  4) .  The  bending  of  the 
recompression  wave  between  items  2  and  3  is  due  to  the 
interactions  between  the  recompression  wave  and  the  initial 
nozzle  expansion  reflecting  off  the  lower  solid  boundary. 
Unlike  Figures  4.3  and  4.4,  the  recompression  in  Figures  4.7 
and  4.8  does  not  reflect  off  the  nozzle  wall.  Due  to  the 
increased  velocity  in  the  nozzle  for  the  Mach  15  flight 
condition,  the  recompression  is  carried  downstream  without 
reflecting  on  the  nozzle  wall. 

Figures  4 . 9  and  4 . 10  show  the  pressure  and  temperature 
distributions,  respectively,  along  the  nozzle  wall  for  the 
FDS  and  SCHNOZ  frozen  flow  solutions  for  the  20.6  deg  nozzle 
attachment  angle  at  the  Mach  15  flight  condition.  Figure 
4.9  shows  that  the  initial  expansion  is  much  more  rapid  for 
the  Mach  15  flight  condition  than  the  Mach  10  flight 
condition  shown  in  Figure  4.5.  This  occurs  because  of  the 
increased  velocity  of  the  Mach  15  combustor  exit  results. 

The  FDS  frozen  flow  solution  calculates  a  lower  pressure  and 
temperature  rise  at  the  recompression.  Again,  the  FDS 
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solution  is  artificially  dampened  by  the  additional 
computational  steps.  Table  4.2  shows  that  the  FOS  solution, 
with  the  grid  packing,  uses  over  twice  as  many  computational 
steps  as  the  SCHNOZ  solution.  Aside  from  the  slight 
differences  at  the  recompression,  the  FDS  and  SCHNOZ 
property  distributions  at  the  nozzle  wall  are  nearly 
identical  for  the  Mach  15  flight  condition  and  20.6  deg 
nozzle  geometry. 

The  static  pressure  contour  plot  for  the  FDS  imperfect 
gas  solution  with  the  20.6  deg  nozzle  attachment  angle  is 
illustrated  in  Figure  4.11  for  the  Mach  20  flight  condition. 
The  static  pressure  contour  plot  for  the  SCHNOZ,  inviscid 
frozen  flow  solution  for  the  same  flight  condition  and 
nozzle  geometry  is  illustrated  in  Figure  4.12.  Figures  4.11 
and  4 . 12  show  a  good  agreement  between  the  two  solution 
methods  for  the  general  characterization  of  the  nozzle 
flowfield.  The  initial  nozzle  expansion  (item  1)  is  more 
severe  for  the  Mach  20  flight  condition  because  of  the 
greater  nozzle  inlet  velocity.  The  recompression  due  the 
change  in  nozzle  curvature  (item  2)  is  quickly  carried 
downstream  and  coalesces  into  a  shock  wave  within  the  nozzle 
( item  3 ) . 

The  calculated  pressure  distribution  along  the  nozzle 
wall  for  both  the  FDS  and  SCHNOZ  frozen  flow  solutions  is 
presented  in  Figure  4.13  for  the  Mach  20  flight  condition 
and  the  20.6  deg  nozzle  attachment  angle.  Figure  4.13  shows 
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that  the  expansion  is  more  rapid  for  the  Mach  20  flight 
condition  than  for  the  Mach  15  flight  condition  (Figure 
4.9).  The  rapid  expansion  is  enhanced  by  the  increased 
nozzle  inlet  velocity  for  the  Mach  20  flight  condition.  The 
FDS  frozen  flow  solution  calculates  a  slightly  lower 
pressure  rise  from  the  recompression  than  the  SCHNOZ 
solution  due  to  the  artificial  damping  induced  by  the  FDS 
grid  packing. 

The  calculated  pressure  distribution  along  the  nozzle 
wall  for  both  the  FDS  and  SCHNOZ  frozen  flow  solutions  is 
presented  in  Figure  4.14  for  the  Mach  10  flight  condition 
and  the  38  deg  nozzle  attachment  angle.  Comparing  Figure 
4.14  to  Figure  4.5,  the  pressure  distribution  plot  for  the 
20.6  deg  nozzle  geometry,  shows  that  the  38  deg  attachment 
angle  causes  a  more  rapid  and  severe  expansion  about  the 
nozzle  attachment  radius.  The  initial  expansion  about  the 
38  deg  attachment  radius  achieves  a  pressure  of 
approximately  one-third  of  that  for  the  expansion  about  the 
20.6  deg  attachment  radius.  The  increased  attachmnet  angle 
allows  for  a  greater  nozzle  expansion.  Figure  4.14  shows 
nearly  a  10%  difference  in  the  FDS  and  SCHNOZ  calculated 
magnitudes  of  the  recompression  pressure  rise.  Table  4.2 
shows  that  the  FDS  solution,  with  grid  packing,  uses  over 
twice  as  many  computational  steps  as  the  SCHNOZ  solution. 

The  artificial  damping  from  the  additional  computational 
steps  is  more  pronounced  for  the  38  deg  attachment  angle 
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nozzle,  as  compared  to  the  20.6  deg  nozzle  geometry,  because 
the  recompression  wave  is  stronger  due  to  a  greater  change 
in  nozzle  curvature  downstream  of  the  nozzle  attachment 
point. 

Figures  4.3  through  4.14  show  a  generally  good 
agreement  between  the  FDS  and  SCHNOZ  frozen  flow  solutions 
in  capturing  the  physics  of  the  hypersonic  propulsive  nozzle 
flowfield.  The  maximum  differnece  in  the  pressure 
distribution  along  the  nozzle  wall  is  approximately  10%, 
occurring  at  the  Mach  10  flight  condition  for  the  38  deg 
nozzle  attachment  angle.  The  greatest  difference  between 
the  two  solutions,  however,  is  in  the  required  computational 
time.  Table  4.2  shows  that  the  FDS  imperfect  gas  method 
requires  CPU  times  between  6.4  and  9.9  times  greater  than 
the  SCHNOZ  frozen  flow  solution. 

4.4  Viscous  Effects 

The  implementation  of  the  viscosity  terms  in  the  SCHNOZ 
code  proved  to  be  ineffective  in  truly  capturing  any  viscous 
effects.  Figures  4.15  and  4.16  show  the  static  pressure  and 
temperature  distributions,  respectively,  along  the  nozzle 
wall  for  the  viscous  and  inviscid  calculations  of  the  38  deg 
parabolic  nozzle  contour  at  the  Mach  10  flight  condition. 
There  is  no  discernable  difference  between  the  property 
distributions  for  the  inviscid  and  viscous  flow 
calculations.  This  same  pattern  was  seen  in  all  other 

nozzle  geometries  and  flight  conditions. 
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Table  4.3  shows  the  difference  in  the  CPU  time  and  the 
calculated  nozzle  thrust  for  the  inviscid  and  viscous  SCHNOZ 
results.  The  differences  in  calculated  nozzle  wall  thrusts 
between  the  viscous  and  inviscid  SCHNOZ  solutions  are  less 
than  0.1%  for  all  flight  conditions  and  geometries. 

Compared  to  the  inviscid  solution,  the  viscous  solution 
required  a  factoral  increase  of  up  to  1.25  times  the 
computational  time.  Essentially,  the  implementation  of  the 
viscous  terms  in  the  SCHNOZ  code  provided  no  additional 
information  to  the  flowfield  solution  at  the  expense  of 
additional  CPU  time.  The  ineffectiveness  of  the  SCHNOZ  code 
to  model  the  viscous  effects  is  attributable  to  the  slip 
flow  boundary  condition  employed  to  remove  the  subsonic 
portion  of  the  boundary  layer. 

4.5  Finite^Rate  Chemistry  Effects 

Frozen  and  finite-rate  thermodynamic  models  were 
considered  at  each  flight  condition  to  study  the  effect  of 
chemistry  on  the  calculated  nozzle  flowfield  for  both  the 
20.6  deg  and  38  nozzle  contours.  Figures  4.17  and  4.18  show 
the  pressure  and  temperature  distributions,  respectively, 
along  the  nozzle  wall  for  the  38  deg  parabolic  nozzle 
contour  for  the  Mach  10  flight  condition.  Figure  4.17  shows 
that  there  is  a  slight  difference  in  the  pressure 
distribution,  as  the  frozen  flow  calculates  a  smaller 
pressure  rise  downstream  of  the  initial  expansion  compared 
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to  the  finite-rate  chemistry  results.  The  difference  in  the 
pressure  distribution  results  in  a  difference  in  the 
calculated  nozzle  thrust.  Table  4.4  shows  the  increase  in 
CPU  time  and  the  percent  change  in  the  calculated  thrust  for 
the  finite-rate  chemistry  runs  compared  to  the  frozen  flow 
results.  The  calculated  nozzle  wall  thrust  is  greater  for 
the  finite-rate  chemistry  SCHNOZ  solutions  than  the  frozen 
flow  solution.  The  increase  in  the  calculated  nozzle  wall 
thrust  is  due  to  the  increase  in  the  pressure  distribution 
along  the  nozzle  wall  for  the  finite-rate  chemistry 
solutions.  The  percent  change  in  nozzle  thrust  is  a  maximum 
of  4%  at  the  Mach  20  flight  condition  and  is  independent  of 
the  two  nozzle  geometries  considered.  The  solution  of  the 
viscous  PNS  equations  including  finite-rate  chemistry 
effects,  required  CPU  run  times  approximately  2.5  times 
greater  than  the  inviscid  SCHNOZ  frozen  flow  computations. 

The  difference  in  the  temperature  distribution  is  more 
pronounced  than  the  difference  in  pressure  distribution 
between  the  two  thermodynamic  models.  Figure  4.18  shows  the 
increased  temperature  distribution  along  the  nozzle  wall  of 
the  finite-rate  chemistry  solution  compared  to  the  frozen 
flow  solution  for  the  Mach  10  flight  condition  and  the  38 
degree  parabolic  nozzle  contour.  The  finite-rate  chemistry 
model  allows  for  dissociation  and  chemical  reactions 
throughout  the  nozzle  flowfield.  As  shown  in  Table  4.5,  the 
mole  fraction  of  the  combustion  products,  H2O,  increases 
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from  the  nozzle  inlet  to  exit.  The  increase  in  H2O 
indicates  that  the  combustion  process  is  continuing  as  the 
flow  expands  through  the  nozzle.  The  combustion  process 
releases  energy  into  the  flowfield.  This  additional  energy 
results  in  a  higher  temperature  for  the  finite-rate 
chemistry  solution  than  the  frozen  flow  solution. 

4.6  FDS  Perfect  Gas  Trends 

4.6.1  FDS  Perfect  Gas  vs .  Imperfect  Gas 
Figures  4.19  and  4.20  present  the  static  pressure  and 
temperature  distributions,  respectively,  along  the  nozzle 
wall  for  the  FDS  imperfect  and  perfect  gas  solutions  to  the 
38  deg  parabolic  nozzle  contour  at  the  Mach  10  flight 
condition.  The  perfect  gas  calculates  both  a  higher 
pressure  and  temperature  distribution  along  the  nozzle  wall. 
The  calculation  of  the  pressure  using  the  linearized- 
approximate  FDS  method  was  virtually  the  same  for  both  gas 
models  considered,  see  Appendix  C.  However,  there  is  a 
greater  difference  in  the  calculation  of  temperature  because 
of  the  variation  of  the  ratio  of  specific  heats,  y,  allowed 
in  the  imperfect  gas  model.  The  variation  of  y  with 
temperature  for  the  imperfect  gas  and  the  increased  modes  of 
molecular  energy  allow  the  imperfect  gas  to  more  equally 
distribute  the  energy  conversion  across  the  recompression 
wave.  The  perfect  gas  absorbs  a  greater  percentage  of  the 
energy  change  into  the  translational  mode  of  energy,  thus 
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resulting  in  an  increase  in  temperature. 

The  addition  of  the  imperfect  gas  model  to  the  FDS 
algorithm  causes  a  significant  increase  in  computational 
time.  Table  4.6  shows  the  increase  in  CPU  time  and  the 
percent  change  in  the  calculated  thrust  for  the  imperfect 
gas  FDS  results  compared  to  the  perfect  gas  FDS  results. 

CPU  run  times  were  up  to  20  times  greater  for  the  imperfect 
gas  model  than  the  perfect  gas  FDS.  The  added  run  time  is 
caused  by  the  additional  iterations  required  in  the  solution 
of  the  Riemann  problem  downstream  of  a  discontinuity  and  the 
decoding  of  the  Riemann  flux  vector  components. 

4.6,2  FDS  Perfect  Gas  vs.  SCHNOZ  Finite-Rate  Chem 

It  is  interesting  to  note  that  the  FDS  perfect  gas 
assumption  results  more  closely  follow  the  finite-rate 
chemistry  results  of  the  SCHNOZ  code  than  do  frozen  flow 
results.  Figures  4.17  and  4.18  show  that  the  SCHNOZ  finite- 
rate  chemistry  model  solution  calculates  slightly  higher 
pressure  and  temperature  distributions  than  the  SCHNOZ 
frozen  flow  solution.  The  higher  pressure  and  temperature 
distributions  occur  because  of  the  added  energy  to  the  flow 
from  the  combustion  process  continuing  throughout  the 
nozzle.  Similarly,  Figures  4.19  and  4.20  show  that  the 
perfect  gas  FDS  solution  calculates  higher  pressure  and 
temperature  distributions  than  the  imperfect  gas  FDS.  The 
higher  pressure  and  temperature  distributions  occur  because 
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of  the  increase  in  the  energy  absorbed  by  the  translational 
mode  of  energy  downstream  of  the  flowfield  discontinuities. 

An  adhoc  comparison  of  the  FDS  perfect  gas  solutions 
and  SCHNOZ  finite-rate  chemistry  model  solutions  is 
presented  in  Figures  4.21  through  4.24.  Figure  4.21 
illustrates  the  pressure  distribution  along  the  nozzle  wall 
for  the  SCHNOZ  finite-rate  chemistry  model  solution  and  the 
FDS  perfect  gas  solution  for  the  20.6  deg  nozzle  geometry  at 
the  Mach  10  flight  condition.  Figure  4.21  shows  a  good 
agreement  between  the  two  calculated  pressure  distributions 
along  the  nozzle  wall.  Figure  4.22  shows  the  temperature 
distribution  along  the  nozzle  wall  for  the  same  flight 
condition  and  geometry.  Except  for  the  differences  at  the 
recompression  peaks,  the  FDS  perfect  gas  temperature 
distribution  is  offset  by  a  constant  amount  from  the 
temperature  distribution  determined  from  the  SCHNOZ  finite- 
rate  chemistry  model. 

Figure  4.23  shows  a  relatively  constant  offset  in  the 
FDS  perfect  gas  temperature  distribution  and  the  SCHNOZ 
finite-rate  chemistry  results  for  the  Mach  15  flight 
condition  and  the  20.6  deg  nozzle  geometr^'.  This  same 
constant  offset  is  evident  in  the  two  calculated  temperature 
distributions  for  the  20.6  deg  nozzle  geometry  at  the  Mach 
20  flight  condition  illustrated  in  Figure  24.  The 
artificial  damping  because  of  the  grid  packing  used  in  the 
FDS  perfect  gas  solution  contributes  to  difference  in  the 
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temperature  distributions  at  the  recompression  peaks  for 
Figures  4.22  through  4.24. 

Table  4.7  presents  the  CPU  time  comparison  for  the  FDS 
perfect  gas  and  SCHNOZ  finite-rate  chemistry  solutions.  The 
perfect  gas  FDS  algorithm  runs  6  to  8  times  faster  than  the 
SCHNOZ  finite-rate  chemistry  runs.  The  constant  offset  in 
the  temperature  distributions  in  Figures  4.22-4.24  indicates 
that  a  slight  change  in  the  y  value  used  for  the  FDS  perfect 
gas  solution  could  yield  a  more  accurate  temperature 
distribution.  The  choice  of  a  perfect  gas  y  that  is 
representative  of  the  actual  flow  physics  would  allow  the 
perfect  gas  FDS  solution  to  provide  valuable  information  on 
the  flow  phenomena  trends  exhibited  in  the  nozzle  with  a 
tremendous  savings  in  computational  time. 
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Table  4.1.  Comparison  of  grid  packing  efficiencies,  FDS 
perfect  gas  solutions,  M=15,  38  deg  nozzle. 


Packing  fraction 

Axial  Steps 

CPU  time  (sec) 

None 

387 

34.1 

1.10 

806 

74.2 

1.05 

1291 

119.6 

Table  4.2.  Computational  planes  and  CPU  time  (seconds) 
comparison  for  the  SCHNOZ  and  FDS  inviscid 
frozen  flow  solutions. 


Flight 

condition 

Axial 

Steps 

FDS 

Axial 

Steps 

SCHNOZ 

CPU 

Time 

FDS 

CPU 

Time 

SCHNOZ 

M=10, 

20.6  nozzle 

1461 

562 

2214.8 

224.8 

M=15, 

20.6  nozzle 

934 

443 

1767.1 

187.2 

M=20, 

20.6  nozzle 

691 

403 

1111.6 

173.2 

M=10, 

38  nozzle 

1258 

561 

2105.2 

224.3 

M=15, 

38  nozzle 

788 

480 

1467.1 

200.4 

Table  4.3.  CPU  time  and  calculated  thrust  comparisons  of 
SCHNOZ  viscous  vs.  inviscid  solutions,  frozen 
flow,  38  deg  parabolic  nozzle  contour. 


Flight  condition 

factoral  increase 
in  CPU  time  vs. 
inviscid  solution 

percent  change  in 
nozzle  thrust  vs. 
inviscid  solution 

Mach  10 

1.14 

0.00 

Mach  15 

1.18 

0.08 

Mach  20 

1.25 

0.02 
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Table  4.4.  CPU  time  and  calculated  thrust  comparisons  of 
SCHNOZ  finite-rate  chemistry  vs.  frozen 
solutions,  38  deg  parabolic  nozzle  contour. 


Flight  condition 

factoral  increase 
in  CPU  time  vs. 
frozen  flow 

percent  increase 
in  nozzle  thrust 
vs.  frozen  flow 

Mach  10 

2.65 

3.89 

Mach  15 

2.47 

3.33 

Mach  20 

2.49 

4.06 

Table  4.5.  Combustion  products,  H2O,  mole  fraction  at 

nozzle  inlet  and  exit  for  SCHNOZ  finite-rate 
chemistry  solutions. 


Flight  condition 

inlet  H2O  mole 
fractions 

exit  H2O  mole 
fractions 

M=10,  20.6  nozzle 

0.2442 

0.2808 

M=15,  20.6  nozzle 

0.2075 

0.2436 

M=20,  20.6  nozzle 

0.1608 

0.1982 

Table  4.6.  CPU  time  and  calculated  thrust  comparisons  of 
FDS  imperfect  gas  vs.  perfect  gas  solutions, 
20.6  deg  parabolic  nozzle  contour. 


Flight  condition 

factoral  increase 
in  CPU  time  vs. 
perfect  gas 

percent  change  in 
nozzle  thrust  vs. 
perfect  gas 

Mach  10 

15.79 

-6.70 

Mach  15 

19.94 

-1.40 

Mach  20 

17.53 

-1.40 
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Table  4.7.  CPU  time  (seconds)  comparison  for  the  FDS 

perfect  gas  and  SCHNOZ  finite-rate  chemistry 
solutions. 


Flight  condition 

FDS 

perfect 

gas 

SCHNOZ 
finite 
rate  chem 

factoral 
increase 
in  time 
vs .  FDS 

M=10, 

20.6  nozzle 

140.2 

702.4 

5.01 

M=15, 

20.6  nozzle 

88.6 

553.2 

6.24 

M=20, 

20.6  nozzle 

63.4 

496.8 

CO 

• 

M=10, 

38  nozzle 

121.6 

679.0 

5.58 

M=15, 

38  nozzle 

75.1 

583.7 

7.77 
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Figure  4.1.  static  pressure  distribution  along  nozzle  wall 
for  different  grid  packing,  FDS  solution,  M=15, 
perfect  gas,  38  deg  parabolic  nozzle. 


Figure  4.2.  static  pressure  distribution  along  nozzle  wall 
for  different  grid  packing,  FDS  solution,  M=15, 
frozen  flow,  38  deg  parabolic  nozzle. 
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Figure  4.3.  Static  pressure  contours  (atm),  FDS  frozen  flow 
solution,  M=10,  20.6  deg  parabolic  nozzle 
contour . 


Figure  4.4.  Static  pressure  contours  (atm) ,  SCHNOZ  inviscid 
frozen  flow  solution,  M=10,  20.6  deg  parabolic 
nozzle  contour. 
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Figure  4.5.  Static  pressure  distribution  along  nozzle  wall, 
FDS  and  SCHNOZ  inviscid  solutions,  frozen  flow, 
M=10,  20.6  deg  parabolic  nozzle  contour. 


x/t> 


Figure  4.6.  static  temperature  distribution  along  nozzle 

wall,  FDS  and  SCHNOZ  inviscid  solutions,  frozen 
flow,  M=10,  20.6  deg  parabolic  nozzle  contour. 
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Figure  4.7.  static  pressure  contours  (atm),  FDS  frozen  flow 
solution,  M=15,  20.6  deg  parabolic  nozzle 
contour . 


Figure  4.8.  static  pressure  contours  (atm),  SCHNOZ  inviscid 
frozen  flow  solution,  M=15,  20.6  deg  parabolic 
nozzle  contour. 
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Figure  4.9.  Static  pressure  distribution  along  nozzle  wall, 
FDS  and  SCHNOZ  inviscid  solutions,  frozen  flow, 
M=15,  20.6  deg  parabolic  nozzle  contour. 


Figure  4.10.  Static  temperature  distribution  along  nozzle 

wall,  FDS  and  SCHNOZ  inviscid  solutions,  frozen 
flow,  M=15,  20.6  deg  parabolic  nozzle  contour. 


Figure  4.11,  Static  pressure  contours  (atm),  FDS  frozen 

flow  solution,  M=20,  20.6  deg  parabolic  nozzle 
contour . 


Figure  4.12.  Static  pressure  contours  (atm),  SCHNOZ 

inviscid  frozen  flow  solution,  M=20,  20.6  deg 
parabolic  nozzle  contour. 
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Figure  4.13.  Static  pressure  distribution  along  nozzle  wall, 
FDS  and  SCHNOZ  inviscid  solutions,  frozen  flow, 
M=20,  20.6  deg  parabolic  nozzle  contour. 


Figure  4.14.  Static  pressure  distribution  along  nozzle  wall, 
FDS  and  SCHNOZ  inviscid  solutions,  frozen  flow, 
M=10,  38  deg  parabolic  nozzle  contour. 
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Figure  4.15.  Static  pressure  distribution  along  nozzle 

wall,  SCHNOZ  viscous  and  inviscid  solutions, 
frozen  flow,  M=10,  38  deg  parabolic  nozzle. 
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Figure  4.16.  static  temperature  distribution  along  nozzle 
wall,  SCHNOZ  viscous  and  inviscid  solutions, 
frozen  flow,  M=10,  38  deg  parabolic  nozzle. 
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Figure  4,17.  Static  pressure  distribution  along  nozzle 

wall,  SCHNOZ  frozen  and  finite  rate  chemistry 
viscous  solutions,  M=10,  38  deg  nozzle. 


Figure  4.18.  Static  temperature  distribution  along  nozzle 
wall,  SCHNOZ  frozen  and  finite  rate  chemistry 
viscous  solutions,  M=10,  38  deg  nozzle. 
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figure  4.19.  Static  pressure  distribution  along  nozzle 

wall,  FDS  perfect  and  imperfect  gas  solutions, 
M=10,  38  deg  parabolic  nozzle  contour. 


4000 


Figure  4.20.  Static  temperature  distribution  along  nozzle 

wall,  FDS  perfect  and  imperfect  gas  solutions, 
M=15,  38  deg  parabolic  nozzle  contour. 
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Figure  4.21.  Static  pressure  distribution  along  nozzle 
wall,  SCHNOZ  finite-rate  chemistry  and  FDS 
perfect  gas  solutions,  M=10,  20.6  deg  nozzle. 


Figure  4.22.  Static  temperature  distribution  along  nozzle 
wall,  SCHNOZ  finite-rate  chemistry  and  FDS 
perfect  gas  solutions,  M=10,  20.6  deg  nozzle. 
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Figure  4.23,  Static  temperature  distribution  along  nozzle 
wall,  SCHNOZ  finite-rate  chemistry  and  FDS 
perfect  gas  solutions,  M=15,  20.6  deg  nozzle. 


Figure  4.24.  Static  temperature  distribution  along  nozzle 
wall,  SCHNOZ  finite-rate  chemistry  and  FDS 
perfect  gas  solutions,  M=20,  20.6  deg  nozzle. 
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V.  Conclusions /Recommendations 


5 . 1  Conclusions 

The  overall  comparison  of  the  two  computational 
methodologies  revealed  that  currently  the  SCHNOZ  code  is 
more  efficient  than  the  FDS  frozen  flow  code  for  analyzing  a 
given  nozzle  inlet  profile  and  geometry.  The  SCHNOZ  and  FDS 
code  demonstrate  a  good  agreement  in  modeling  the  flow 
physics  in  the  hypersonic  nozzle.  The  inclusion  of 
viscosity,  turbulence,  and  finite  rate  chemistry  models 
should  make  the  SCHNOZ  code  more  accurate.  However,  the 
implementation  of  the  solid  wall  slip  boundary  condition  to 
remove  the  subsonic  portion  of  the  viscous  boundary  layer 
fails  to  truly  model  the  effects  of  viscosity  in  the  nozzle 
flowfield.  The  chemistry  effects  in  the  nozzle  flowfield 
have  a  small  (less  than  5%)  effect  on  the  overall  nozzle 
thrust.  However,  even  these  small  thrust  differences  are 
significant  in  the  optimized  design  of  a  NASP-type  vehicle 
and  the  inclusion  of  finite  rate  chemistry  effects  should  be 
incorporated  into  the  computational  solution. 

The  FDS  code  utilizing  the  perfect  gas  assumption  is 
useful  and  efficient  in  optimizing  a  nozzle  or  cowl  geometry 
for  a  given  inlet  condition.  The  perfect  gas  FDS  results 
show  the  proper  trends  in  flow  phenomena  and  flowfield 
properties,  and  provide  a  good  ballpark  design  to  analyze  by 
a  more  accurate  method.  The  added  CPU  time  of  the  more 
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accurate  FDS  frozen  flow,  compared  to  the  perfect  gas,  is 
attributable  to  the  additional  iterations  required  in  the 
solution  of  the  Riemann  problem  and  the  flux  vector  decode. 
This  added  CPU  time  greatly  reduces  the  efficiency  of  the 
FDS  method.  An  improved  iteration  scheme  for  the  frozen  gas 
thermodynamic  model  for  the  FDS  could  should  markedly 
improve  its  performance  in  terms  of  CPU  time. 

5 . 2  Recommendations 

According  to  the  SCHNOZ  operator's  manual  (20) ,  a  new 
version  of  the  SCHNOZ  code  was  to  be  developed  that  includes 
a  viscous  sublayer  model  and  an  option  for  equilibrium 
chemistry.  Future  computations  of  the  nozzle  flowfield 
should  be  run  with  this  version  of  the  SCHNOZ  code  or  a 
different  PNS  code  that  includes  the  viscous  sublayer  model 
to  better  capture  the  viscous  flow  effects.  With  proper 
consideration  of  the  viscous  effects  on  the  solution,  a 
study  to  determine  whether  the  PNS  equations  are  necessary 
to  model  the  flow  physics  would  prove  beneficial. 

Continual  improvements  in  computational  methodologies 
and  computational  capabilities  need  to  be  investigated  and 
analyzed  to  meet  the  demanding  challenges  of  the  hypersonic 
propulsive  nozzle  environment.  The  current  version  of  the 
SCHNOZ  program  is  not  optimized  for  use  on  a  supercomputer, 
and  the  FDS  code  is  not  vector izable  because  of  the  nature 
of  the  solution  to  the  Riemann  problem.  Future  studies 
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should  be  made  into  the  efficiency,  in  terms  of 
computational  time  and  cost,  for  a  vector izable  PNS  code  run 
on  a  supercomputer. 

The  true  accuracy  of  the  computer  codes  can  only  be 
determined  by  a  comparison  to  exact  solutions.  In  lieu  of 
unlikely  actual  flight  test  data,  an  unclassified 
experimental  analysis  of  a  simple  nozzle  geometry  should  be 
undertaken  to  validate  the  findings  and  trends  of 
computational  runs.  Also,  computational  runs  using  actual 
combustor  exit  data  would  prove  to  be  of  merit  in  comparison 
to  previously  assumed  or  derived  nozzle  inlet  profiles. 
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Appendix  A:  Governing  Equations 


A.l  Parabolized  Navier-Stokes  (PNS)  Equations 

The  set  of  PNS  equations  implemented  in  the  SCHNOZ  code 
were  derived  from  the  full  Navier-Stokes  equations  using 
time-averaging  of  density,  pressure,  and  viscous  stresses 
and  mass-averaging  of  the  remaining  variables  (15:2-1). 

Using  Reynold's  averaging  and  ignoring  unsteady  terms  and 
streamwise  derivatives  of  viscous  stress  terms,  the 
turbulent  Navier-Stokes  equations  were  parabolized  (4:907). 
Ignoring  turbulent  stress  terms  in  the  normal  momentum 
equation  and  diffusion  terms  in  the  energy  equation,  because 
of  an  assumed  unity  Lewis  number  (4:907),  the  resulting 
planar  PNS  equations  are  obtained: 

Continuity : 

d[pu)  ^  d{pv)  ^  Q 
dx  dy 

X-Momentum: 

d[p  +  pu^)  ^  d(puv)  ^  d  f-  du\ 
dx  dy  3yV  dy) 

Y-Momentum: 

d(puv)  ^  d{p+pv^)  ^  d  hr  dv\ 
dx  dy  3yV  ^7/ 


(A.l) 


(A. 2) 


(A. 3) 


A-1 


Total  Enthalpy: 


djpuH)  ^  djpvH) 
dx  dy 


_d_ 

dy 


f— 1 

Pr  dy 

[  2  / 

d  I  ^  dH' 
dy\  Pr  dy 


(A.  4) 


Full  details  on  the  derivation  of  these  equations  are 
available  in  the  work  of  Sisilian  (14)  and  their  application 
to  the  SCHNOZ  code  is  presented  by  Dash  (20) . 

For  chemically  reacting  viscous  flows  a  finite-rate 
chemistry  model  is  used,  and  an  additional  equation  is  added 
to  the  PNS  equation  set  to  account  for  the  various  species 
continuity; 


.  a{pvai)  _  d(  p  da^] 

d  dy  “i  ■  dy) 


(A. 5) 


where  the  subscript  i  represents  the  individual  species 
considered. 


A. 2  Turbulence  Model 

Due  to  computer  limitations  and  the  small  space  scales 
of  turbulent  motion,  the  effects  of  turbulent  disspation  and 
diffusion  are  only  computationally  feasible  using  turbulence 
models  (9:659).  The  SCHNOZ  program  contains  two  turbulence 
models,  the  kc  and  kW  high  Reynolds  number  versions.  The 
turbulence  model  used  in  this  study  is  the  common  two- 
equation  ke  eddy  viscosity  version  which  models  the 
turbulent  kinetic  energy,  k,  and  energy  dissipation,  6.  The 
turbulent  viscosity,  is  determined  from: 
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(A.6) 


l^t 


where  the  length  scale  parameter  used  is  the  turbulence 
energy  dissipation  rate,  e  (5:507).  The  energy  dissipation 
rate  is  related  to  k  and  the  turbulent  length  scale,  £, 
through  the  following  relation  (15:10): 


k2/2 


(A.7) 


The  parabolized  form  of  the  turbulence  modelling  equations 
are  developed  from  the  full  Navier-Stokes  equations  for  the 
production,  transport,  and  dissipation  of  k  and  e  with  the 
omission  of  all  axial  derivatives  (5:507).  For  planar  flow 
the  turbulence  equations  are  then  given  by: 


and  : 


_e 

k 


(A. 8) 


(A.  9) 


The  constants  used  in  the  parabolized  turbulence  modelling 
equations  are  (15:10): 

C,  =  1.43  Uk  =  1.0 

Cj  =  1.92  a,  =  1.3 

=  0.09 

The  ke  model  is  developed  from  incompressible 
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assumptions  and  overestimates  mixing  rates  for  high  Mach 
number  flows,  such  as  those  in  a  hypersonic  nozzle  (5:507). 
Therefore,  a  compressibility  correction,  K(M^)  ,  is  employed 
in  the  SCHNOZ  code  to  arrive  at  a  corrected  value  of  the 
turbulent  viscosity  given  by: 

(A.  10) 

where  is  the  characteristic  Mach  number  of  the  turbulence 
(5:507).  This  characteristic  Mach  number  of  the  turbulence 
is  found  from: 

Af-  =  (A. 11) 

^max 

where  is  the  maximum  value  of  k  at  each  axial  location 
and  a  ^  is  the  local  speed  of  sound  at  each  k^  (5:507). 

The  values  of  K(M^)  are  determined  from  comparison  with 
experimental  data  and  are  shown  in  Figure  A.l  (5:507). 

A. 3  Euler  Equations 

The  Euler  equations  are  the  zero  viscosity  and  zero 
conductivity  limits  of  the  Navier-Stokes  equations  (9:659). 
The  effects  of  viscosity  and  the  conduction  of  heat  are  not 
included  in  the  Euler  equations.  The  Euler  equations  as 
utilized  in  the  FDS  computer  program  of  Doty  are: 
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Continuity : 


(A. 12) 


(A. 13) 


(A.14) 


(A. 15) 


Traditionally,  solutions  of  the  Euler  equations  require 
significantly  less  computational  time  than  comparable  PNS 
equation  solutions  (2:236). 
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Appendix  B:  Thermodynamic  Model 


B . 1  Perfect  Gas 

The  thermodynamic  model  employed  greatly  influences  the 
determination  of  flow  properties  within  the  hypersonic 
nozzle.  The  complexity  of  the  thermodynamic  models  depend 
on  the  degree  to  which  the  model  actually  captures  the  true 
physics  of  the  flow  (13:51).  The  simplest  model  is  the 
perfect  gas.  A  calorically  perfect  gas  model  assumes  no 
intermolecular  forces  and  constant  specific  heats,  c^  and  Cp, 
resulting  in  a  constant  specific  heat  ratio,  7.  The  perfect 
gas  obeys  the  thermal  equation  of  state  (1:381): 

p  =  pRT  (B.l) 

where  R  is  the  gas  constant  determined  by  the  molecular 
weight,  MW,  of  the  gas  and  the  universal  gas  constant, 

R  =  (b.2) 

MW 

A  perfect  gas  assumes  no  change  in  the  molecular  composition 
of  the  gas.  Hence,  R  remains  constant  since  there  is  no 
mechanism  to  change  the  MW  of  the  gas.  The  enthalpy  and 
internal  energy  of  a  calorically  perfect  gas  are  functions 
only  of  temperature  (1:388). 
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(B.3) 


h  =  h{Tl  =  CpT 

Q  =  Q(T)  =  C^T  (B.4) 

The  total  internal  energy  for  a  calorically  perfect  gas  is 
given  by: 

pe  =  (B.5) 

The  calorically  perfect  gas  assumption  is  the  most 
widely  used  assumption  in  flow  analyses.  However,  the 
calorically  perfect  gas  assumption  breaks  down  when  high 
temperatures  are  encountered  within  the  flowfield.  The 
perfect  gas  does  not  account  for  the  internal  structure  of 
the  molecule,  nor  does  it  allow  for  the  dissociation  of 
gaseous  molecules  with  increasing  temperature.  These 
shortcomings  are  adressed  by  the  calorically  imperfect  gas 
model  and  the  finite-rate  chemistry  model,  respectively. 

B.2  Calorically  Imperfect  Gas 

The  next  step  in  complexity  is  the  calorically 
imperfect  gas.  The  calorically  imperfect  gas  still  assumes 
that  intermolecular  forces  are  negligible  and  obeys  the 
thermal  equation  of  state.  Eg.  (B.l).  The  basic  difference 
between  the  calorically  imperfect  gas  and  perfect  gas  is  in 
the  recognition  of  the  proper  internal  structure  of  the 
gaseous  molecule  and  the  modes  of  internal  energy  which  may 
be  excited.  The  differences  in  the  two  assumptions  is 
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illustrated  by  considering  air  as  a  diatomic  gas  and 
ignoring  the  relatively  small  contributions  of  the 
electronic  mode  of  energy  (18:136). 

In  general  the  specific  internal  energy,  now  denoted  by 
e  and  not  ti,  for  a  gas  is  additive  and  represented  by 
(18:128)  : 


^total  ^tzans  *  .  \i  ^intern  (  •  ) 

incern 

where  e;^^  represents  the  internal  modes  of  energy 
available.  For  a  diatomic  gas  such  as  air,  the 
translational,  rotational,  and  vibrational  modes  of  energy 
are  available,  even  at  relatively  moderate  temperatures 
(18:222).  Using  a  statistical  mechanics  derivation  Vincenti 
and  Kruger  provide  the  following  expressions  for  the 
translational,  rotational,  and  vibrational  internal  modes  of 
energy  for  a  diatomic  gas  (18:133-135): 


^tzana  2 

(B.7) 

^zot  ~ 

(B.8) 

(B.9) 

where  0^  is  the  characteristic  temperature  for  vibration. 

The  perfect  gas  assumption  allows  only  for  the  translational 
and  rotational  internal  modes  of  energy.  The  specific 
internal  energy  for  a  perfect  gas  is  then  given  by: 
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^perfect  ^trana  ®rot  2^^  (B.IO) 

The  calorically  imperfect  gas,  accounting  for  the  internal 
energy  mode  of  vibration  is: 

^ total  ~  ^ttana  *  ®rot  ^  ®vii>  ~  "2  ^  ®vii>  ( B .  1 1 ) 

Comparing  Eqs.  (B.IO)  and  (B.ll)  shows  that  the  energy  of 
gaseous  molecules  is  spread  over  more  modes  of  energy  for  a 
calorically  imperfect  gas  than  for  the  perfect  gas.  This 
difference  influences  the  determination  of  flow  properties 
downstream  of  discontinuities.  For  perfect  gas  flow  across 
a  shock  wave,  the  kinetic  energy  of  the  flow  in  front  of  the 
shock  is  converted  to  translational  and  rotational  energy 
downstream  of  the  shock.  For  the  calorically  imperfect  gas, 
some  of  this  kinetic  energy  conversion  may  be  absorbed  by 
the  vibrational  mode  of  energy,  thus  decreasing  the  amount 
of  energy  absorbed  by  the  rotaional  and  translational  modes. 
Since  temperature  is  a  direct  measurement  of  kinetic 
(translational)  energy,  the  increased  energy  absorbed  in  the 
translational  mode  of  the  perfect  gas,  compared  to  the 
calorically  imperfect  gas,  causes  the  perfect  gas  model  to 
overpredict  the  temperature  downstream  of  the  shock  (1:510). 
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B.3  Complex  Chemical  Mixtures 

The  combustion  products  of  a  hydrogen-fueled  air- 
breathing  SCRAM jet  engine  are  a  complex  mixture  of  gaseous 
species  consisting  largely  of  water  and  nitrogen.  The 
specific  heat  for  a  mixture  of  n  gaseous  species  is  given  by 
(1:387)  : 

Cp  -  <B-“) 

where  Cj  is  the  mass  fraction  of  the  species.  The  specific 
heat  of  the  species,  Cpj,  for  a  mixture  of  thermally  perfect 
gases  is  dependent  upon  the  temperature  of  the  mixture 
(1:388).  Using  curve  fitted  JANAF  thermochemical  data,  the 
least  squares  coefficients  of  specific  heats,  enthalpy,  and 
entropy  parameters  may  be  obtained  as  functions  of 
temperature.  In  terms  of  the  least  squares  coefficients,  Cp 
is  written  as  (13:56): 

Cp^  =  (a  +  bT  +  cT^  +  dT^  +  er^)i2^ggCj  (B.13) 

Enthalpy  for  a  thermally  perfect  gas  mixture  is  given 
by  (1:397): 

h  =  *  [CpdT  {B.14) 

*0 

Substituting  Eg.  (B.13)  into  Eg.  (B.14)  yields  the  least 
sguares  coefficients  form  of  enthalpy  for  species  i  (13:57): 

-  [{/2o  +  aT  +  ^  £t^  ir  ^  *  ^T^)Rgas  ^ 
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where  hT.29g  is  the  enthalpy  of  the  species  at  the  reference 
temperature  of  298  K. 

Entropy  for  a  thermally  perfect  gas  is  given  by 
(1:396) : 


s 


dT 


-  Rln^ 
Po 


=  <b  -  2?ln-^ 
Po 


(B.16) 


Again,  substituting  Eg.  (B.13)  into  the  integral  portion  of 
Eg.  (B.16)  yields  the  least  sguares  coefficient  form  of  the 
species  entropy  parameter,  <pi,  for  a  species  i  (13:57): 

4)^  =  (4)0  +  ar  +  +  j)Rgas<^i  (b.i?) 

The  total  specific  heat,  enthalpy,  and  entropy  of  the 
thermally  perfect  mixture  is  then  found  from  the  summation 
of  the  individual  species  value.  As  implemented  in  the 
thermally  perfect  gas  model  (calorically  imperfect) ,  the 
chemical  composition  of  the  mixture  remains  constant,  or 
frozen,  throughout  the  nozzle. 


B.4  Finite-Rate  Chemically  Reacting  Flow 

The  finite-rate  chemistry  kinetics  model  accounts  for  a 
reacting  mixture  of  gaseous  species  and  the  dissociation  of 
gaseous  molecules  at  high  temperatures.  Dissociation  alters 
the  molecular  composition  and  molecular  weight  of  the 
gaseous  mixture  and  occurs  around  2000  K  for  Oj,  and  begins 
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at  4000  K  for  Nj  (1:451).  The  change  in  molecular 
composition  and  molecular  weight  alters  the  properties  of 
the  flow  and  the  determination  of  flowfield  properties 
downstream  of  discontinuities. 

For  a  chemically  reacting  mixture  of  n  different 
gaseous  species,  Xj,  the  general  chemical  reaction  equation 
is: 


(B.18) 


where  p/  and  I'j”  are  the  stochiometric  coefficients  of  the 
reactants  and  products,  respectively.  The  stochiometric 
coefficients  are  positive  for  products  and  negative  for 
reactants.  The  constants  kf  and  are  the  forward  and 
reverse  reaction  rate  constants.  The  net  rate  of  production 
of  a  given  species  i  is  given  by  (1:493): 


d[X^] 

dt 


vf 


(B.19) 


where  [X,]  denotes  the  concentration  of  species  i  in  moles 
per  unit  volume.  The  rate  constants,  kf  and  k^,  are  related 
by  the  equilbriura  constant  based  on  concentration,  K^,  such 
that  (1:493) : 


(B.20) 


where  the  concentration  based  equilibrium  constant,  K^,  can 
be  found  from  the  equilibrium  constant  based  on  partial 
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pressures,  Kp,  using  the  equation: 


\  ■^uoll 


t) 


(B.21) 


The  equilibrium  constant  based  on  partial  pressures,  Kp,  for 
a  given  chemical  reaction  is  determined  from  the  Gibbs  free 
energy  and  temperature  (1:406): 


K  (  T)  *  G  Rfjj^yT 


(B.22) 


In  Eg.  (B.22)  AC***  is  the  Gibbs  free  energy  of  the  products 
minus  the  Gibbs  free  energy  of  the  reactants  for  a  given 
chemical  reaction  with  all  species  evaluated  at  a  pressure 
of  one  atmosphere  and  a  specified  reference  temperature 
(1:402-406).  As  implemented  in  the  SCHNOZ  computer  program, 
values  of  are  determined  from  a  thermodynamic  datafile 
in  the  computer  code. 

The  species  concentration,  [X;] ,  in  Eg.  (B.19)  can  also 
be  written  as 


[^iJ  =  pF, 


(B.23) 


where  p  is  the  density  of  the  mixture  and  Fj  is  the  mass 
fraction  of  species  i  divided  by  its  molecular  weight: 


F,  = 


(B.24) 


The  net  production  rate  equation.  Eg.  (B.19),  can  then  be 
written  as 
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The  chemical  reactions  implemented  in  the  finite-rate 
chemistry  model  of  the  SCHNOZ  code  for  this  study  are: 

H  +  O2  OH  +  O 

H2  +  O  OH  +  H 

OH  +  OH  H2O  +  O 

H2  *  OH  ^  H2O  *  O  ^5^26 

O  +H+M'>*‘  OH  +  M 
H  +  H  +  H2  +Af 

H  +  OH  +  H2O  +  M 

O  +0+Af-*02  +  M 


For  the  reactions  considered  in  Eq.  (B.26)  M  is  a  third  body 
which  can  be  a  molecule  of  any  species  present  in  the 
mixture.  The  finite-rate  chemistry  model  in  the  SCHNOZ  code 
assumes  that  all  third  bodies  have  equal  efficiencies  in 
contributing  to  the  reactions  (15:2-9).  The  rate  of 
production  of  a  species  i  for  the  system  is  then  the  sum  of 
the  rates  of  production  of  species  i  for  the  individual 
reactions  considered  (1:496-497). 
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Appendix  C:  Flux-Difference-Splitting 


C . 1  Introduction 

The  flux-difference-splitting  (FDS)  method  is  a 
technique  capable  of  capturing  complicated  flows  with  strong 
property  gradients.  The  FDS  method  requires  the  solution  of 
the  Riemann  problem  at  a  given  axial  location  to  advance 
downstream  in  a  marching  fashion.  FDS  takes  advantage  of 
the  wave-like  nature  of  the  Riemann  problem  to  split  the 
flow  vector  fluxes  along  proffered  paths  of  propagation 
(6:153).  Full  details  on  the  solution  of  the  Riemann 
problem  for  planar  supersonic  flow  assuming  a  thermally  and 
calorically  perfect  gas  are  provided  in  reference  (6) .  This 
solution  procedure  was  modified  for  a  calorically  imperfect 
gas  assumption  in  reference  (13).  The  FDS  method  is 
independent  of  the  Riemann  problem  solution  method. 

C.2  The  Riemann  Problem 

The  Riemann  problem  is  represented  in  Figure  C.l. 
Godunov  proposed  that  the  general  flow  property  ,  St', 
distribution  can  be  modelled  as  a  series  of  uniform  flow 
regions  with  a  discontinuity  occuring  half-way  between  the 
nodes  of  interest  (6:11).  In  Figure  C.l  the  solid  line  ' 
represents  the  arbitrary  ’i'  distribution  and  the  dashed  line 
represents  the  Godunov  regions  of  uniform  flow. 
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The  collapse  of  the  discontinuity  to  the  midpoint, 
j+1/2,  results  in  waves  that  are  propagated  along 
characteristics  based  on  the  wave  angle  (6:9)  as  depicted  in 
Figure  C.2.  Waves  (1)  and  (3)  can  be  any  combination  of 
expansion  or  compression  waves,  and  wave  (2)  is  a  contact 
surface  that  separates  regions  (2)  and  (4) .  Flow  property 
discontinuities  are  present  across  the  waves,  however  the 
contact  surface,  wave  (2) ,  cannot  support  a  pressure  nor 
flow  angle  discontinuity  (6:163).  The  values  of  flow 
properties  in  regions  (0)  and  (6)  are  known  from  the  initial 
values  or  the  solution  of  the  previous  Riemann  problem. 

C.3  Solution  of  the  Riemann  Problem 

The  solution  of  the  Riemann  problem  is  the 
determination  of  the  primitive  flow  variables,  p,  p,  u,  and 
V,  in  regions  (2)  and  (4)  as  shown  in  Figure  C.2.  Doty  (6) 
details  three  different  methods  that  may  be  used  to  solve 
the  Riemann  problem.  The  first  is  an  exact  solution 
procedure  where  iterations  of  non-linear,  coupled  equations 
are  required  for  a  non-isentropic  compression  wave  and 
iterations  of  the  non-linear  Prandtl-Meyer  relations  are 
necessary  for  simple  expansions.  Further  iterations  are 
required  to  match  the  pressure  in  regions  (2)  and  (4)  to 
satisfy  the  contact  surface  boundary.  The  second  method  is 
an  exact-approximate  solution.  This  is  similar  to  the  exact 
method  but  treats  compression  waves  as  isentropic.  The 
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isentropic  Prandtl-Meyer  relations  may  then  be  used  to  solve 
for  flow  variables  across  any  wave.  The  non-linear  Prandtl- 
Meyer  relations  still  require  iterative  solution  techniques 
as  well  as  further  iterations  of  the  entire  solution  to  meet 
the  contact  surface  pressure  requirements.  The  third  method 
is  the  linearized-approximate  solution  and  is  used 
exclusively  in  this  study. 

The  linearized-approximate  solution  treats  all  waves  as 
isentropic  expansions  and  compressions  which  may  be  solved 
by  the  Prandtl-Meyer  relations.  Furthermore,  the  Prandtl- 
Meyer  relations  are  linearized  to  allow  for  a  closed  form 
algeabraic  equation  requiring  no  iterations.  The 
linearized-approximate  solution  to  the  Riemann  problem  is 
dependent  on  the  assumed  thermodynamic  model.  The 
determination  of  the  primitive  variables  in  region  (2)  and 
(4)  is  more  complicated  for  the  calorically  imperfect  gas 
and  requires  an  iterative  technique  for  the  solution  of 
temperature  . 


C.3.1  Linearized-Approximate  Solution,  Perfect  Gas 
For  isentropic,  planar,  steady  flow,  the  differential 
form  of  the  compatibilty  relations,  valid  along  Mach  lines 
is  (6:165) : 

v/Af^-1  dp  ±  pv^dd  =  0  (c.i) 

where  the  (+)  sign  refers  to  positive  characteristics  and 
the  (-)  sign  refers  to  negative  characteristics.  Using  the 
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definitions  of  velocity  magnitude,  V  =  +  v^,  and  flow 

angle,  6  =  atan(v/u) ,  and  the  following  relationship,  p  = 
7P/a^,  Eq.  (C.l)  can  be  rearranged  to  yield: 

±  ±Uill^d(v/u)  .  0  (C. 

P  yjM^-1 

Introducing  the  following  definitions; 

^  _  (yu^/a^) 

yfM^ 

a  =  v/u 

d[ln{p)]  = 

P 

a  more  efficient  form  of  the  compatibility  equations  may  be 
written  as: 

d[ln(p)]  ±  (z)do  =  0  (C.6) 

Linearizing  Eq.  (C.6)  yields: 

A[ln(p)]  ±  {z)Ao  =  0  (C.7) 

Referring  to  Figure  C.2,  a  positive  wave,  wave  (3),  is 
required  to  pass  information  from  region  (0)  to  region  (2) . 
Using  the  (+)  sign  for  the  postive  wave  from  Eq.  (C.7)  and 
using  the  regions  (0)  and  (2)  in  the  difference  operator 
gives: 

([ln{p)]2-  [ln{p)]o)  +  (Zq)  (Oj-Oq)  =0  (C.8) 


(C.3) 

{C.4) 

(C.5) 
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Rearranging  Eg.  (C.8)  yields: 


[ln(p)]2  +  (Zo)®2  =  [ln(p)]o  +{Zo)®o  (C.9) 

Similarly,  a  negative  wave  is  required  to  pass  information 
from  region  (6)  to  (4)  and  from  Eg.  (C.7)  is  expressed  by: 

[ln(p)]4  -  {^5)04  =  [ln(p)]6  (C.IO) 

The  pressure  and  flow  angle  (or  slope,  a)  are  required  to 
match  across  the  contact  surface.  Thus: 

O4  =  O2  (C.ll) 


P4  =  P2  (C.12) 

Equations  (C.9),  (C.IO),  (C.ll),  and  (C.12)  now  constitute  a 
set  of  four  equations  with  four  unknowns.  After 
substitution  and  rearranging  the  four  equations,  an 
expression  for  the  solution  of  the  flow  angle  in  region  (4) , 
04,  is  obtained  in  terms  of  known  properties  in  regions  (0) 
and  (6) : 


[In  (p)]o-[ln(p)e-»-(z6)P6'*'(^o)°o 
(Ze+^o) 


(C.13) 


Equation  (C.IO)  is  then  solved  for  P4: 

P4  =  exp(  [ln{p)  +  26(04-06))  (C.14) 
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The  density  and  speed  of  sound  in  regions  (2)  and  (4)  are 
determined  from  isentropic  ratios  across  waves  (1)  and  (3) : 


Po 


(C.15) 


^2  =  [YPa/Pzf^^ 


(C.16) 


^  = 
Pe 


(C.17) 


^4  =  [YP4/p4f^^  (C.18) 

Conservation  of  stagnation  enthalpy  across  waves  (1) 
and  (3)  is  used  to  obtain  the  velocity  components  in  regions 
(2)  and  (4) .  Recall  from  Appendix  B  that  for  a  calorically 
perfect  gas  h  =  CpT.  Thus,  across  wave  (3)  : 


K . 
■ 

Equation  (C.19)  can  be  solved  explicitly  for  M4;  and  the 
velocity  components,  in  terms  of  the  flow  slope,  are 
(13:62) : 

U4  =  ^434003  (tan‘^  (04)  ]  (C.20) 

V4  =  W4a4sin(tan'^  (04)  ]  (C.21) 

A  similar  procedure  is  used  to  find  the  primitive  flow 
variables  in  region  (2) . 
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C.3.2  Linearized-Approximate  Solution,  Imperfect  Gas 
The  determination  of  the  pressure  in  region  (4) ,  Eq. 
(C.14)  remains  valid  for  a  calorically  imperfect  gas.  To 
determine  the  remaining  primitive  flow  variables  a  relation 
is  needed  that  does  not  rely  on  a  calorically  perfect  form 
of  the  conservation  of  stagnation  enthalpy. 

As  outlined  in  reference  (13) ,  the  determination  of 
temperature  uses  relations  that  take  advantage  of  the 
assumed  isentropic  nature  of  the  flow.  For  a  thermally 
perfect  gas,  the  differential  quantity  of  entropy  in  a 
system  is  (1:397); 

ds  =  Cp^  -  R-^  (C.22) 

Integration  yields; 

s  =  =  <t>  -  -R  In-S.  (C.23) 

i  ^  T  Po  Po 

where  the  entropy  parameter,  <p,  is  defined  as  (13:63): 

*  •  (C.24) 

Substituting  Eq.  (B.13)  from  Appendix  B  for  the  specific 
heat,  Cp,  for  a  mixture  of  gases,  into  Eq.  (C.24)  yields: 

<|>  =  |<|>o  +  aInT  +  JbT  +  +  -^jR  (C.25) 

A  change  in  entropy  across  wave  (3)  is  found  from  Eq.  (C.23) 
and  written  in  terms  of  the  entropy  parameters  in  regions 
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(4)  and  (6)  as: 

S4  -  Se  =  i|>4  - 

For  isentropic  flow,  Eq.  (C.26)  reduces  to 

4),  .  4 

where  <p^  is  found  using  Eq.  (C.25)  based  on  the  known 
parameters  in  region  (6).  Once  Eq.  (C.27)  is  solved  for  <p^, 
Eq.  (C.25)  may  be  iteratively  solved  for  the  temperature,  t. 
The  thermal  equation  of  state  is  then  used  to  solve  for  the 
density  in  region  (4) . 

The  conservation  of  stagnation  enthalpy  is  used  in  a 
different  form  than  for  the  perfect  gas  to  determine  the 
velocity  components  in  region  (4) .  Conservation  of 
stagnation  enthalpy  across  wave  (3)  provides  for  the 
determination  of  the  stagnation  enthalpy  in  region  (4) : 

(C.28) 

The  static  enthalpy,  h^,  is  a  function  of  Tj  and  is  found 
using  the  methods  described  in  Appendix  B  and  the  velocity 
components,  \x^  and  v^,  are  known.  Using  the  definition  of 
the  slope,  a  =  v/u,  Eq.  (C.28)  can  be  rewritten  as: 
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(C.29) 


ht,  =  K  *  -|u|(l+o|) 

Solving  for  the  axial  velocity  component  in  region  (4) ,  U4 
gives: 


U4  = 


N 


I+O4 


(C.30) 


The  radial  component  of  velocity,  V4,  is  then  computed  using 
the  slope  in  region  (4) : 


V4  =  U4O4 

The  Mach  number  is  then  computed  : 

V 

M.  = 


(C.31) 


(C.32) 


C.4  Riemann  Fluxes 

After  solving  the  Riemann  problem,  the  primitive 
variables  are  combined  to  form  the  Riemann  flux  vectors  in 
regions  (0) , (2) ,  (4),  and  (6)  (6:171).  The  E  and  F  flux 
vectors  are  presented  in  section  3 . 1  and  are  repeated  here 
for  convenience. 


pu 

pv 

pu^+p 

F  = 

pVTi 

puv 

u(pe+p). 

/  * 

pv2+p 

v(pe+p). 

(C.33) 


The  Riemann  fluxes  are  calculated  for  each  of  the  components 
of  the  E  and  P  vectors.  For  example,  the  first  component  of 
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the  E  vector,  El  is 

pu 

and  is 

computed 

in  each  of  the 

Riemann  regions  (0) , 

(2) 

,  (4), 

and  (6) : 

iEDo 

=  Po^o 

(C. 

34) 

=  P2U2 

(C. 

.35) 

{El)^ 

=  P4U4 

(C. 

,36) 

(El)e 

=  Pe^e 

(C. 

,37) 

The  remaining  Riemann  flux  vector  components  are  found  in  a 
similar  manner. 


C.5  Calculation  of  the  Flux  Differences 

The  Riemann  flux  differences  are  the  differences  in  the 
Riemann  flux  vector  components  taken  across  waves  (1),  (2) 
and  (3) .  For  example,  the  difference  of  the  Riemann  fluxes 
for  the  first  component  of  the  E  vector,  dEl,  across  waves 
(1)  ,  (2)  and  (3)  are: 


{EDo 

(C.38) 

=  (SI),  - 

{ED  2 

(C.39) 

=  (Bl),  - 

{El)i 

(C.40) 

The  sum  of  the  total  contributions  across  all  three  waves 
gives  the  total  contribution  at  the  midpoint  of  the  Riemann 
region,  node  j+1/2  (6:182): 
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The  remaining  Riemann  flux  vector  components  are  computed  in 
a  similar  manner. 

C.6  Splitting  the  Flux  Differences 

Referring  to  Figure  C.3  the  information  is  known  at 
node  j  on  axial  plane  i,  and  the  solution  is  sought 
downstream  at  node  j  on  axial  plane  i+1.  The  solution  of 
the  Riemann  problem  calculated  the  fluxes  at  the  midpoints 
j+1/2  and  j~l/2.  The  Riemann  flux  differences  were 
calculated  and  are  now  split  along  their  direction  of 
propagation  to  determine  which  information  from  Riemann 
nodes  j+1/2  and  j-1/2  will  reach  node  j  at  plane  i+1. 

Following  the  method  detailed  by  Doty  (3:183-186)  the 
flux  differences  are  propagated  in  a  direction  determined  by 
streamlines  and  characteristics.  At  node  j+1/2,  if  the 
slope  of  waves  (1),  (2),  or  (3)  is  negative,  the  Riemann 
flux  differences  across  those  waves  contribute  to  the 
solution  at  the  downstream  node  (i+l,j)  (6:185). 

Conversely,  at  node  j+1/2,  if  the  slope  of  waves  (1),  (2), 
or  (3)  is  positive,  the  Riemann  flux  differences  across 
those  waves  do  not  contribute  to  the  solution  at  (i+l,j) 
but,  would  contribute  to  the  solution  at  node  (i+l,j+l). 
Similarly,  at  Riemann  node  j-1/2,  if  the  slope  of  waves  (1) , 
(2) ,  or  (3)  is  positives.  The  Riemann  flux  differences 
across  those  waves  contribute  to  the  solution  at  node 
(i+l,j).  The  Riemann  flux  differences  are  essentially  split 
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into  positive  and  negative  contributions  at  each  Riemann 
node. 

The  contributions  at  the  downstream  node  (i+l,j)  are 
found  by  summing  the  positive  contributions  of  the  Riemann 
flux  differences  from  node  j-1/2  and  the  negative 
contributions  of  the  Riemann  flux  differences  from  node 
j+1/2; 

dB/.y,  =  [d£7“i7/l*  *  *  Kj”.7fl*  «=•«> 

where  the  positive  or  negative  signs  denote  positively  or 
negatively  sloped  flux  differences.  It  is  possible  that 
all,  some,  or  none  of  the  Riemann  flux  differences  at  a 
Riemann  midpoint  node  may  contribute  to  the  downstream 
solution  at  a  Riemann  node,  depending  on  their 
characteristic  slope.  The  same  procedure  is  applied  to  the 
F  vector.  Further  details  on  the  splitting  of  the  flux 
differences  are  presented  in  Appendix  J  of  reference  (6) . 
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Figure  C.l.  General  property  distribution  and  Riemann 
description  (6:135). 


Figure  C.2.  Riemann  problem  for  planar  supersonic  flow 
and  resulting  wave  pattern  (6:174). 
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Appendix  D:  FDS  Decode  Procedures 


D . 1  Introduction 

The  upwind  FDS  algorithm  solves  the  downstream  E  and  P 
flux  vectors  based  on  the  solution  of  the  Riemann  problem 
generating  Riemann  fluxes  that  are  then  split  along  their 
physically  correct  path  of  propagation  as  discussed  in 
Appendix  C.  To  continually  march  the  solution  from  plane  i 
to  i+1  the  values  of  the  primitive  variables  are  required  at 
each  plane  to  solve  the  Riemann  problem.  A  decode  procedure 
is  required  at  each  axial  plane  to  extract  the  primitive 
variables  from  the  calculated  flux  vectors.  This  decode 
procedure  is  dependent  upon  the  thermodynamic  model  utilized 
in  the  flowfield  solution. 

D.2  Perfect  Gas  Decoding 

The  perfect  gas  decode  procedure  is  detailed  in 
reference  (6)  and  summarized  herein.  The  four  components 
of  the  E  vector,  El,  E2,  E3,  and  E4,  are  known  from  the 
spatial  marching  of  the  FDS  algorithm.  The  decoding  of  the 
E  vector,  see  Eg.  (3.2),  requires  the  simultaneous  solution 
of  five  equations  for  the  five  unknowns  p,  u,  v,  p,  and  e 
contained  in  the  E  flux  vector. 

The  pe  term  contained  in  the  E4  component  can  be 
written  in  terms  of  the  internal  energy  and  flow  energy; 
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(D.l) 


pe  =  pQ  +  -|p  (u^+v^) 

A 

For  a  perfect  gas  the  specific  internal  energy,  a,  is 
written  as  (1:388): 


a  =  c^T 


(D.2) 


where  c,  for  a  perfect  gas  can  be  written  in  terms  of  the 
gas  constant,  R,  and  the  ratio  of  specif c  heats,  7,  as: 


(D.3) 


Using  the  thermal  equation  of  state  and  the  relations  for  ti 
in  Eg.  (D.2)  and  c^  in  Eg.  (D.3), the  expression  for  the  term 
pe  in  Eg.  (D.l)  can  be  rewritten  as: 


pe  =  +  .|p{u2+^2) 


(D.4) 


Using  Eg.  (D.4)  for  the  pe  term,  the  E4  component  is 


given  by: 


E4  =  u[(-^  +  -|p(u2  +  v2))  +  p] 
Y  1  ^ 


(D.5) 


Expanding  Eg.  (D.5)  and  rearranging  yields; 


E4  =  upf— +  -ipu(u2+v2 


(D.6) 


The  E2,  E3,  and  E4  components  can  be  written  in  terms  of  the 
El  component,  pu,  to  give: 


E2  =  {El)  u  *  p 


(D.7) 
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E3  =  {El)v 


(D.8) 


E4 


+  1{E1)  (U2)  +  ^{eD  (v2) 


(D.9) 


Solving  Eq.  (D.8)  for  the  radial  velocity,  v,  gives: 


V  = 


E3 

El 


(D.IO) 


Similarly,  Eq.  (D.8)  can  be  solved  for  the  pressure,  p: 

p  =  E2  -  {El)  U  (D.ll) 

Eqs.  (D.IO)  and  (D.ll)  can  now  be  substituted  into  the  E4 
component  equation,  Eq.  (D.9)  to  yield  an  expression  that  is 
solely  a  function  of  the  known  E  vector  components  and  one 
unknown ,  u : 


E4 


[  {E2)  -  {ED  u  +  1  {ED  (u2)  +  -1  {ED 

\E3-\ 

lY-lJ 

2  2 

El 

(P.12) 


Casting  this  expession  as  a  quadratic  equation  yields: 


2(y-1) 


{ED 


- 

{E2)' 

u  + 

L  Y-1  J 

{E4)  - 


1 

2 


{E3)^ 

{ED 


=  0 


(P.13) 


This  quadratic  equation  is  then  used  to  find  the  value  of 
the  axial  velocity,  u: 


_  -b  -  4ac 

2a 


(D.14) 


where : 


a 


2(y-1) 


{ED 


(D.15) 
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(D.16) 


b  =  (E2) 

Y-1 


(D.17) 


The  remaining  primitive  flow  variables  are  then  found  from: 


V  = 


E3 

El 


(D.18) 


p  =  E2~  (El)  U 


(D.19) 


P  = 


El 

U 


(D.20) 


pe  =  ^+-|p(u2+v2) 


(D.21) 


The  perfect  gas  decode  procedure  is  then  complete. 

D.3  Imperfect  Gas  Decoding 

The  imperfect  gas  decode  procedure  development  follows 
the  procedure  presented  in  refernece  (13) .  The  difference 
in  the  decoding  of  the  E  vector  for  the  calorically 
imperfect  gas  is  that  the  internal  energy  term,  pe,  in  the 
E4  component  must  be  consistent  with  the  imperfect  gas 
assumption  (13:75).  Specific  internal  energy,  Ci,  for  the 
calorically  imperfect  gas  is  given  by  (1:395): 

Q  =  h  -  RT  (D.22) 
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The  total  internal  energy  for  the  imperfect  gas  is  then 
found  by  substituting  Eg.  (D.22)  into  the  pe  expression,  Eg. 
(D.l)  to  yield: 


pe  =  ph  -  p  +  p-|  (u2+v2) 

Substituting  Eg.  (D.23)  into  the  E4  component  gives: 


E4  =  puh  +  ~pu{u^+v^) 

A 


(D.23) 


(D.24) 


Substituting  the  El  component,  pu,  and  the  E3  component, 
puv,  into  Eg.  (D.24)  gives  the  following  form  of  the  E4 


component : 


E4  =  h{El)  +  -|{^l)u2  + 


(D.25) 


Eg.  (D.25)  is  now  a  function  of  the  axial  velocity,  u,  and 
the  enthalpy,  h.  The  enthalpy  is  a  function  of  temperature 
which  is  determined  from  the  least  sguares  coefficient  as 
discussed  in  Appendix  B.  The  thermal  eguation  of  state, 
which  is  still  valid  for  the  calorically  perfect  gas,  and 
the  El  component  are  combined  to  yield  an  expression  for  u: 


y  ^  (El)  RT 


(D.26) 


Solving  the  E2  component  for  p  and  substituting  the  result 
into  Eg.  (D.26)  yields: 


(El)  RT 
E2-(E1)  u 


(D.27) 
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Rearranging  Eg.  (D.27)  a  quadratic  for  u  is  developed: 

iEl)u^  +  {E2)U  +  (E1)RT=  0  (D.28) 

which  when  solved  for  u  yields: 

=  E2  +  y/{E2)  -4  (El)  (  {EIYRTT  (D.29) 

2  (El) 

The  E4  component  is  now  a  function  only  of  temperature  with 
both  h  and  u  dependent  on  the  temperature  value.  The 
temperature,  T,  is  found  iteratively  using  the  secant  method 
and  Eqs.  (D.28)  and  (D.29)  (13:77). 

Once  the  temperature  is  determined,  u  is  determined 
from  the  quadratic  relation.  The  remaining  primitive 
variables  are  found  using  the  perfect  gas  decode  procedure 
of  Eqs.  (D.18)  -  (D.21) . 
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